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1.0 


INTRODUCTION 


This  document  Is  a  final  report  on  studies  of  high  current  beam 
propagation  In  low  density  gases.  The  study  has  extended  over  several  years 
and  has  resulted  In  two  Informal  Interim  topical  reports  —  resummarized 
herein.  Three  topics  have  been  addressed:  Ohmic  Hall  currents,  non-local 
non-Ohmlc  conductivity  modelling,  and  non-linear  electromagnetic  algorithms. 
The  Hall  and  non-local  conductivity  studies  are  Important  primarily  for  high 
beam  currents  In  low  density  gases;  the  non-linear  electromagnetic  field 
solvers  are  Important  for  studying  beams  In  low  density  channels. 

Our  Initial  studies  of  Hall  current  effects  were  motivated  by  the 
observation  that  the  EfxB  drift  of  plasma  electrons  would  tend  to  expand  the 
radial  distribution  of  plasma  current  and  reduce  the  normal  axial  current. 

The  Initial  supposition  was  that  these  stabilizing  effects  would  be  most 
Important  for  high  currents  and  low  channel  densities,  since  the  Hall  currents 
scale  directly  with  the  Larmor  frequency,  r.^,  and  Inversely  with  the  collision 
frequency,  y  .  Countervailing  effects,  however.  Include: 

o  The  high  level  of  current  neutralization  characteristic  of  high- 
current  beams, 

o  Saturation  of  Hall-current  effects  by  their  associated  potential 
buildup,  and 

o  Transition  In  the  high-current  regime  to  Spltzer  conduction  In  the 
beam  core,  enhancing  on-axis  conductivity  relative  to  neutral- 
col  1  Islon-domlnated  conductivity  In  the  wings  of  the  beam's 
profile. 

Further  complexities  In  the  actual  situation  Include  the  Interplay 
between  E  -  and  Ez~f1elds  In  the  progress  of  avalanche  breakdown  at  the  head 
of  the  beam,  the  effects  of  air  chemistry  processes  (especially  recombina¬ 
tion),  and  the  role  of  radiative  cooling  In  determining  the  onset  of  Spltzer 
conduction.  Progress  In  modelling  these  effects  Included  the  development  of  a 
dependable  electromagnetic  algorithm;  the  algorithm  Included  forward-time 
differencing  and  an  Iterative  approach  to  find  fields  In  the  highly  nonlinear 
situation  created  by  the  presence  of  tensor  conductivity  effects. 


Studies  of  the  effects  of  Hall  currents  In  a  local.  Ohmic  approxi¬ 
mation  are  reported  In  Section  2.  No  consequences  significant  to  beam  propa¬ 
gation  or  stability  are  seen  for  gas  pressures  greater  than  76  torr.  At  lower 
pressures,  however,  the  usual  local.  Ohmic  air  conductivity  models  employed  In 
the  stuay  are  Inappropriate.  Significant  effects  which  must  be  considered 
Include: 


o  Time-dependent  (delayed)  and  nonlocal  Ionization  production.  In 
contrast  to  the  usual  local.  Instantaneous  production  models, 

o  Contributions  to  the  total  current  from  the  delta-ray  components 
of  the  beam-driven  Ionization  cascade, 

o  Non-Ohmlc  conductivity  effects, 

o  Transport  of  electrons  from  their  local  volume  of  creation,  and 

o  Non-Maxwel 1  Ian  energy  distribution  of  the  conduction  electrons, 
notably  In  the  "runaway"  regime;  low  effective  collision 
frequency.  In  comparison  to  .  ,  could  especially  enhance  the 
contributions  of  a  high  energy  conduction  group. 

A  conductivity  model  Incorporating  the  above-noted  effects  has 
been  developed  and  Is  reported  In  Section  3.  The  model  treats  three  electron- 
energy  groups:  a  relativistic  (delta)  group,  and  high-  and  low-energy  con¬ 
duction  groups.  The  beam,  as  well  as  each  group,  populate  lower-energy  groups 
via  Ionization  energy  loss;  E-fleld  acceleration  can  move  low-energy  electrons 
upward  to  the  Intermediate  group.  Fluid  models  represent  the  cynamlcs  of  the 
two  lower-energy  groups,  and  a  detailed  air-chemistry  reaction  scheme  also 
modifies  the  population  of  the  low-energy  group.  The  air  chemistry  was  cali¬ 
brated  against  the  detailed  HICHEM  code  In  an  applicable  regime. 

Calculations  with  the  phenomenological  model  sketched  above  showed 
considerably  enhanced  and  broadened  effective  current  profiles  at  times 
between  1  and  10  nsec  for  high-current  beams.  All  of  the  features  of  Improve¬ 
ment  appeared  to  make  significant  contribution  to  the  results  obtained. 

In  the  most  recent  portions  of  our  work,  reported  In  Section  4, 
several  electron agnet 1c  flelc  algorithms  for  beam-propagation  work  wore 
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those  relevant  to  hose  stability)  Into  the  nonlinear  regime  —  1.e.»  to  be 
applicable  for  displacements  comparable  to  beam  radii.  Several  fast  and 
accurate  algorithms  were  found.  All  were  based  on  a  modal  expansion  of  the 
transverse  spatial  dependence  of  the  fields,  and  all  used  an  Iterative  pro¬ 
cedure  to  simplify  the  solution  for  the  mode  amplitudes.  (The  modes  are  all 
coupled  through  the  non-ax  1  symmetry  of  the  conductivity,  and  require.  In 
principle,  simultaneous  solution  for  all  modes.)  Formulations  Involving 
explicit  consideration  of  fields  as  well  as  potentials  In  Lee's  approximation 
(an  important  s  1  rr.pl  If  Icatlon  of  the  frozen-field  approximation)  and  a  frozen- 
field  formulation  were  used.  The  simplified  approaches  were  found  to  agree 
well  with  one  another,  and  with  the  few-polnt  data  available  from  other 
nonlinear  field-solver  algorithms.  Moderate  differences  —  especially  In 
electric  fields  —  were  seen  In  the  frozen-field  formulation  results  for  a 
standard  test  case.  Substantial  Improvements  In  speed  for  the  algorithms 
developed  over  direct-solution  methods  for  the  fully-coupled-mode  equations 
appear  to  have  been  realized. 
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2.0  HALL  CURRENT  EFFECTS  IN  A  LOCAL  CONDUCTIVITY  MODEL 

2.1  Introduction 

This  section  summarizes  our  studies  of  the  effects  of  Ohmic  Hall 
currents  on  the  propagation  of  high  current  beams  In  low  density  gases.  Non- 
Ohmlc,  non-local  Hall  effects  are  discussed  In  Section  3.  Beams  of  very  high 
current  have  attracted  Increased  attention  recently  due  to  potential  stability 
and  application  advantages  (Refs.  1*  2). 

The  Hall  currents  derive  from  the  £x§  drift  of  the  plasma  elec¬ 
trons.  Thus*  due  to  the  pinching  azimuthal  magnetic  field  B^,  plasma  elec¬ 
trons  driven  by  the  Inductive  longitudinal  electronic  field  will  drift  toward 
larger  radii*  and  those  driven  by  the  Coulomblc  radial  electric  field  will 
tend  to  cancel  the  axial  plasma  current.  Initial  estimates  of  the  resulting 
redistribution  of  plasma  currents  relative  to  beam  currents  suggested  poten¬ 
tially  significant  consequences  to  beam  stability  and  propagation.  The 
purpose  of  the  study  herein  reported  was  to  substantiate  the  Initial  estimates 
and  to  explore  some  of  the  consequences. 

Hall  currents  enter  Maxwells  equations  through  terms  proportional 
to  the  Larmor  frequency  divided  by  the  collision  frequency  vm» 

‘=  V-m 

Thus*  Hall  effects  vary  directly  with  net  current  and  inversely  with  gas  den¬ 
sity.  One  anticipates  Hall  effects  to  become  Important  at  very  high  currents 
and  low  densities.  In  the  Ohmic  approximation,  however*  these  two  require¬ 
ments  tend  to  be  mutually  exclusive;  high  current  beams  In  low  density  gases 
tend  to  be  strongly  current  neutralized  and  develop  very  small  net  current  -  ' 
remains  small. 

In  the  body  of  the  pulse,  Hall  effects  achieve  a  quasi-statlonary 
state  and  all  quantities  Important  to  beam  dynamics  become  nearly  Independent 
of  no  matter  how  large  i  may  become.  Thus,  Hall-dependent  consequences 
require  that  current  neutralization  be  substantially  less  than  complete  in  the 
very  early  portions  of  the  pulse. 


For  a  realistic  range  of  high  current  beam  parameters  -  assuming 
Bennett  beam  profiles  and  local  conductivity  models  -  Hall  current  effects 
remain  small.  Near  complete  current  neutralization  obtains  near  the  pulse 
head  and  i  remains  small.  Only  by  artlflcally  constraining  the  conductivity 
profile  evolution  have  we  demonstrated  Important  Hall  current  sensitivities. 
Results  are  thus  very  sensitive  to  non-local  conductivity  effects*  described 
In  Section  3. 

In  Section  2.2  the  general  development  of  Hall  currents  Is  pre¬ 
sented,  Including  monopole  and  dipole  decomposition,  and  Section  2.3  describes 
some  numerically  stable  algorithms  for  monopole  Hall  currents  effects. 

Results  are  In  Section  2.4  and  our  conclusions  to  date  are  In  the  concluding 


2.2 


Theory 


In  this  section  we  derive  the  relevant  electromagnetic  field  equations 
as  employed  in  this  study,  and  discuss  some  important  properties  of  tensor  con¬ 
ductivity.  We  begin  with  Maxwell  equations  in  free  space  (Gaussian  units): 

1  -"-B  =  7  -fI-  +  —  J  :  Ampere's  Law  (2.1) 

C  yl  C 

\>i  =  -  7  :  Faraday's  Law  (2.2) 

C  dt 

where  the  net-current  J  consists  of  both  the  primary  beam-current  and  the 
plasma-current  3  , 

3  =  \  +  3p  .  (2.3) 


For  the  purposes  of  discussion,  a  simple  model  of  the  Hall  current 
effects  may  be  obtained  as  follows.  We  assume  that  the  plasma  electrons  obey 
the  Lorentz  equation  of  motion  with  a  phenomenological  momentum  transfer 
frequency  v 


dv  c*  e  *  „• 
m  -rr=eE  +  -v-B-m  v  v 
e  dt  c  e  m 


(2.4) 


In  tne  local  anprox r  ation,  tr.e  inertial  term  may  be  ignored,  dv/dt  0.  With 
this  approximation,  we  obtain  the  generalized  ohmic  relation  for  the  plasma- 
current  densi ty : 


en  v  =  oE-  —  >-J 
e  v  p 

m  r 

(2.5a) 

e"  n 

e 

m 

e  m 

(2.5b) 

e  B 
ill  c 

(2.5c) 
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o  is  the  usual  scalar  conductivity  and  |H|  is  the  electron  gyrofrequency.  It  is 
the  cross-product  on  the  right  side  of  Eq.  (2.5a)  that  produces  the  Hall  current. 
Later  in  this  section  we  discuss  a  more  rigorous  treatment  where  the  "momentum 
transfer  frequency"  in  Eq.  (2.4)  is  itself  a  function  of  the  plasma  electron 
distribution. 


Solving  for  3p  in  Eq.  (2.5a),  we  find 


3p  =  (1  +  J 
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In  cylindrical  coordinates  (r,4,z),  Eq.  (2.6)  becomes 
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This  form  for  Jp  may  now  be  substituted  back  into  Maxwell  equations  (2.1  and  2.2) 
to  provide  a  closed  set  of  equations  for  generating  (£,£)  from  J.. 


To  make  the  problem  computationally  tractable,  we  decompose  Maxwell 
equations  into  monopole-dipole  subsets.  We  assume  the  following  form  for  E,  £ 
and  j: 


E„  = 


E°  +  E'  e1^ 
r  r 


=  E°  +  E‘  e1^ 


E»  .  E'  e1  ■ 


B  =  B°  +  B '  e1 '  ,  etc. 
r  r  r 


(2.8) 
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where  we  have  made  the  dependence  on  $  explicit.  Furthermore,  we  introduce  the 
retarded  variable  x  =  c  t  -  z  which  we  define  to  measure  the  distance  back  from 
the  pulse  head  (in  the  laboratory  frame  of  reference)  for  a  beam  moving  along  the 
z-axis  with  v2  %  c.  Finally,  we  make  the  frozen  field  assumption;  all  dependent 
variables  are  functions  only  of  the  independent  variables  (x,r),  achieved  by  the 
following  transformation  on  the  partial  derivatives  in  Eqs.  (2.1)  and  (2.2): 


—  c  — 

at  z  ax  z 


(2.9) 


The  monopole  equations  reduce  to  three  equations  relating  E°,  B°,  and 

E°  and  a  second  set  of  three  equations  relating  B°,  E°,  and  B°.  We  may  ignore 

the  second  set  for  jP  =  0  (8°  =  E°  =  B°  =  0).  We  are  left  with  the  set  of 

b$  '  r  0  z 

equations  for  the  axi-symmetric  beam: 
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For  the  dipole  fields,  we  have 
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Note  that  the  T*B  =  0  condition  is  treated  as  an  initial  condition  which  is  pre¬ 
served  by  both  the  monopole  and  dipole  equations. 
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(2. 12d) 


The  dipole  plasma  current-densi ty  also  follows  from  Eq.  (2.7).  After  lineariz¬ 
ing  with  respect  to  the  dipole  variables,  we  find 
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where 


*  =  iTF  ^  Bz> 


(2.13d) 


(2. 13e) 


(2. 13f ) 


This  completes  the  specification  of  the  tensor  plasma  currents  in  terms  of  the 
electromagnetic  field  variables. 

In  the  presence  of  the  Hall  current  terms,  the  radial  plasma  current 

J°  (2.12a)  flows  until  E°  ^  o°E°/on,  in  contrast  to  scalar  calculations  where  ' 

Pr  0  r  6  2 

the  condition  for  J°  ^  0  is  E°  ^  0.  Thus  E°  actually  reverses  in  sign  from  its 

pr  r  r 

original  value  set  up  by  the  Coulomb  field  of  the  beam  particles.  When  Er  reaches 
this  "quasi-equilibrium"  value,  the  expression  for  the  axial  plasma  current 

r  ^ 

as  given  by  Eq.  (2.12c)  becomes 


.o  o  ro 
J  %  o  E 
pz  z 


(2.14) 


which  is  exactly  the  same  as  the  scalar  result  without  Hall  currents.  Similarly, 


-*o  -*o 

for  the  Ohmic  heating  term  J  •  E  ,  we  find 

P 


F-t°  %  o°  E°‘ 


(2.15) 


for  Jpr  %  0,  which  is  the  usual  scalar  relation.  Thus,  the  non-zero  equilibrium 
value  of  E°  in  the  tensor  case  exactly  compensates  for  the  fact  that  o°  <  o°  in 
Eq.  (2.12).  E°  f  0  effectively  describes  a  polarization  of  the  plasma. 
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We  might  conclude  from  this  cancellation  of  terms  that  the  Hall  current 
effects  on  an  axi-symmetric  beam  are  negligible  when  the  radial  plasma  currents 
are  small.  This  would  apply  to  the  body  of  a  beam  pulse  where  charge  neutraliza¬ 
tion  is  complete  and  the  beam  envelope  has  settled  down  to  some  quasi-static 
equilibrium  profile.  However,  if  the  Hall  effects  are  sufficient  to  alter  the 
evolution  of  the  conductivity  channel  and  EM-fields  near  the  beam  head,  where  the 
condition  ^  0  is  not  satisfied,  we  might  still  expect  to  see  a  residue  of  the 
tensor  conductivity  effect  further  back  in  the  beam  where  J°r  %  0.  In  other  words, 
both  o°  and  E°  would  be  noticeably  different  in  both  Eq.  (2.14)  and  (2.15)  in  the 
body  of  the  pulse  due  to  Hall  current  effects  in  the  pulse  head.  It  has  been  the 
purpose  of  the  present  study  to  explore  the  non-equilibrium  evolution  of  tensor 
conductivity  near  the  beam  head. 

The  tensor  conductivity  contribution  to  the  plasma  current  in  Eq.  (2.12) 
relies  on  the  parameter  ft/v  .  Initial  estimates  of  tensor  effects  were  based  on 
the  observation  that  Q/vm  varies  roughly  as  p-1  where  p  is  the  ambient  gas  density 
due  to  the  basic  density  dependence  of  However,  the  gyrofrequency  Q  depends 
linearly  on  the  net  current  which  in  turn  depends  indirectly  (but  strongly)  on  the 
channel  density  and  channel  profile.  Hence,  there  is  no  simple  or  unique  relation¬ 
ship  between  fl/v  and  density  . 

For  a  high  intensity  beam  in  a  low  density  medium,  a  further  complica¬ 
tion  results  from  the  fact  that  the  dominant  contribution  to  the  momentum  transfer 
frequency  vm  shifts  from  electron-neutral  collisions  to  Coulomb  collisions.  Two 
important  changes  occur  at  this  point:  (1)  The  conductivity  becomes  independent 
of  electron  density  (in  fact  decreases  if  rnultiply-charqed  ions  contribute  a  sub¬ 
stantial  fraction  of  the  electrons);  and  (2)  o  scales  as  T3/2  instead  of  roughly 

*1/2  .  ^ 

Tg  .  The  transition  to  the  Coulomb  regime  on-axis  sharpens  the  conductivity 

profile,  allows  more  plasma  current  to  flow  near  the  axis,  and  tends  to  decrease 

the  net  current  flowing  in  this  region.  This  sensitivity  to  electron  temperature 

points  up  the  importance  to  radiative  cooling  effects  --  only  roughly  modeled  here. 


We  conclude  this  section  be  noting  that  the  exact  cancellation  of  the 
tensor  effects  in  Eq.  (2.12)  for  J°r  n,  0  is  a  model  dependent  result.  Its 
validity  requires  that 


(2.16) 


for  all  values  of  B.  A  more  general  theory  in  which  the  random  motion  of  the 
electrons  is  taken  into  account  does  not  satisfy  this  requirement.  For  example, 
if  electron-electron  scattering  is  neglected  (Ref.  3), 


(2.17) 


where 

x  s  j  mv2/kT 


and  v(v)  is  the  electron  collision  frequency  with  neutrals  or  ions.  Only  if  v  is 
independent  of  velocity  v,  or  for  other  very  special  cases,  can  the  relation  (Eq. 
2.16)  be  satisfied  exactly. 
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2.3 


Monopole  Tensor  Electromagnetic  Algorithm 


2.3.1  Introduction 


The  basic  monopole  eguations  in  the  frozen  approximation  (Eg.  2.10)  can 


be  written  in  the  form 


?r  *  4  r  Jz  =  FsT  (r  B)  -  4  "  jb 


(2.18a) 


(2.18b) 


(2.18c) 


where  the  superscript  "o"  has  been  dropped  through  this  entire  section.  [Note 
that  in  the  frozen  approximation,  the  step  Ax  along  the  pulse  length  is  equiva¬ 
lent  to  a  timestep  At  as  the  beam  passes  a  given  fixed  spatial  point.]  The 

plasma  current  densities  J  and  J  take  the  form 

z  r 

Jz  -  °i,  EZ  +  Er 


If  c  HO  (no  Hall  current)  the  equations  are  much  less  strongly  coupled  together, 
and  are  linear  except  for  the  indirect  dependence  of  o  on  t  through  the  ohmic 
heating.  However,  if  Hall  effects  are  included  a  more  explicit  non-linearity 
enters  because  c  is  roughly  proportional  to  B. 


Several  different  numerical  treatments  of  these  equations  have  been 
investigated  in  this  work.  The  algorithm  currently  in  use  is  described  below, 
followed  by  discussions  of  some  other  methods  which  were  tried  and  the  reasons 
they  were  found  to  be  unsati sfactory. 


A  major  difficulty  with  standard  electromagnetic  algorithms  in  the 
charged  particle  beam  context  is  that  the  conductivity  becomes  so  large  that 
it  is  impractical  to  adhere  to  the  timestep  constraint  4-r  c1(At<<  1  which  arises 

in  the  usual  difference  treatment  of  these  equations.  Although  forward  center¬ 
ing  can  eliminate  the  numerical  constraint, i t  can  also  prejudice  the  solution. 
This  is  an  uncertainty  which  has  yet  to  be  completely  resolved  satisfactorily. 

In  at  least  one  case,  however,  our  forward-centered  methods  have  beer,  checked 
by  choosing  the  time  step  to  satisfy  At<<Vr;£;  no  significant  difference  was 
observed. 


2.3.2  The  Present  Algorithm 


Equation  (2.18a)  is  written  in  the  time-integrated  form 


=  Ezo  exp{-4TTj, At)  + 


c  3_ 
r  3r 


( rB )  -  4-rtJg  -  4ttoh 


(l-exp(47ron£t)) 


4tt  a„ 


(2 


where  Ezq  is  the  value  of  Ez  at  the  beginning  of  the  timestep  At,  but  all  other 
quantities  are  evaluated  at  the  end  of  the  timestep.  Although  the  centering  is 
not  perfect  when  4tt  o  At  1,  the  forward  centering  is  found  to  contribute 
significantly  to  stability  and  is  simpler.  When  4tt  aH  »  3/9t,  as  is  usually  the 
case  in  the  spatial  region  occupied  by  the  beam,  the  centering  used  above  gives 
the  correctly-centered  result 


4  TT 
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(c 


E 
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+  0  Er) 


—  -r—  (r  B) . 
r  9r  '  ' 


(2 


In  this  equation  all  quantities  are  evaluated  at  the  forward  time.  It  is  exactly 
the  same  result  as  is  obtained  by  neglecting  displacement  current  in  Eq.  (2.18a), 
and  requiring  that  the  resulting  equation  hold  at  t  +  At. 


Equation  (2.18b)  is  written  as 


from  which 

cr(l  +  4,o, it:  =  (f£  ♦  4,o.  Ez)  it  +  Er0 


(2.22) 


or 


Er  = 


‘ »  ( If  +*=  Ez) 


'ro 


1  +  4' 


:t 


1  +  4;:-.  At 


(2.23) 


In  this  equation,  Er0  is  the  value  at  time  t,  and  all  other  quantities  are  to  be 
evaluated  at  t  +  At.  In  practice,  3B/3t  is  evaluated  as  (B(t+At)  -  B(t))/At,  a 
possible  weakness  to  the  method  because  of  the  centering.  Note  that  the  exponential 
method  as  used  for  the  Ez  equation  could  have  been  used  here  also. 


Equation  (2.23) is  then  substituted  into  Eq.  (2.18c),  written  in  the  form 


o  Ez) 


(2.24) 


and  interpreted  as  being  at  time  t  +  At.  The  result  is 


c;| 

1  +  4-,o|i  At 


(2.25) 


When  this  equation  is  simplified  a  very  important  cancellation  occurs  in  the 
coefficient  of  c  Ez>  which  takes  it  from  being  of  order  unity  to  order  (l+4i  -  At)  , 
a  small  number  when  4iro||At»  1.  An  analogous  cancellation  occurs  if  the  exponential 
form  is  used.  The  resulting  equation  is 
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dr 


4tto„ 


1  +  4i  0||  At 


dB 

dt 


+ 


4t  o  Ez 

1  +  4  7i  o,j  A  t 


(2.26) 


The  final  form  used  for  differencing  is 


3Ez  4- 


(2.27) 


where  again  the  subscript  "o"  indicates  terms  evaluated  at  time  t.  All  other 
terms  are  evaluated  at  t  +  Lt. 

The  two  Equations  (2.19)  and  (2.27)  must  be  solved  simultaneously. 
Since  neither  is  linear  because  of  the  strong  and  direct  dependence  of  oi  on  B, 
an  iterative  procedure  is  used.  The  c  E  /o  term  in  Eq.  (2.27)  is  written  as 

o,  E^  a  B 

°n  cn 

and  a  /0||  is  adjusted  during  the  course  of  the  iteration  process.  Thus  the  co¬ 
efficient  of  B  in  Eq.  (2.27)  becomes 


(i-!M 

\  mc\.  / 


in  the  simple  theory  of  the  Hall  effect.  This  quantity  is  >  1  except  when  Ez  >  0 

at  the  end  of  the  pulse.  In  Eq.  (2.19)  the  term  in  a  E  is  updated  during  the 

course  of  the  iteration  process,  but  is  regarded  as  given  for  each  iteration.  Also 

the  weaker  dependence  of  ;  on  B  (for  most  sit,,ritions  encountered)  is  ignored 

within  each  iteration  step. 

With  these  approximations,  Eqs.  (2.19)  and  (2.27)  can  be  combined  into 
a  single  equation  for  B 


or  r 


^<rB>  -I'  B(1 


(B-Er)0  4 


c  Tf  +  4" v"  TtT  .r 


%-(^.'t(JB+,  Er).  -  E_, 


(2. 28) 
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in  wh  ch 
and 


X  =  exp(-4no  At) 

;■  =  (1  -  x)/(^"  At). 

[Note  that  the  e  Er  term  in  Eq.  (2.19)  could  be  rewritten  as  o  e  Bo  ,  in 
analogy  to  the  treatment  of  the  j  E2  term  in  Eq.  (2.27).  However,  the  explicit 
B-dependence  introduced  in  that  way  would  spoil  the  diagonal-dominant  nature  of 
the  final  equation  for  B.  Also,  E^  itself  has  a  strong  B  dependence  which  would 
be  ignored  by  the  procedure  in  any  case.]  With  E  ,  6  /c  ,  and  c  regarded  as 
fixed,  this  equation  is  solved  for  B  as  a  tridiagonal  difference  system.  The 
values  of  o^E  ,  cjc,  ,  and  o„  are  then  updated  and  the  iteration  procedure  is 
continued  until  the  values  of  the  fields  stabilize.  Although  there  is  no 
guarantee  that  this  iteration  procedure  will  converge,  it  has  been  found 
to  converge  well  in  the  cases  studied. 

This  algorithm  has  been  installed  and  used  in  the  HICHEM  air  chemistry 
code  (Ref.  4)  and  in  the  HIGAP  monopole-envelope  code  (Ref.  5).  It  appears  to 
be  numerically  stable  whether  or  not  o^/o,  is  large.  If  is  set  =  0,  this 
algorithm  (and  minor  variations  of  it)  yields  results  essentially  identical 
(within  reasonable  limits  for  numerical  studies)  to  the  previous  scalar  algorithms 
used  in  HIGAP  and  HICHEM.  The  latter  two  algorithms  are  different  from  each 
other  and  from  the  one  described  above  in  some  aspects  of  the  time  centering. 

The  agreement  of  numerical  results  for  the  scalar  case  does  not  guarantee  the 
tensor  results.  Some  aspect  of  the  centering  or  differencing  method  could  tend 
to  suppress  or  enhance  the  potential  non-linear  evolution  of  the  tensor  solutions 
away  from  the  scalar  solutions. 

2.3.3  Alternative  Numerical  Approaches 

In  the  treatment  discussed  in  Section  2.3.2  above  it  was  necessary  to 
split  up  the  1  and  ||  parts  of  the  plasma  current,  treating  the  JL  part  as  a 
source  term  on  the  right  hand  side  of  the  equation.  An  alternative  treatment 

which  avoids  the  problem  is  to  define 


■ '  s. 


^  Ez  +  1  Er 


./  =  Jz  +  i  Jr. 


then  multiply  Eq.  (2.18b)  by  i  and  add  to  (2.18a).  The  result  is 


Ir*'  •?!?<*>-*■  JB  +  1i 


(2.29) 


Now,/  =  oV,  with  the  definition 


G*  =  0  -  1G 


Then  the  exponential  treatment,  or  full  forward  centering  as  indicated  below, 
can  be  used  to  perform  the  time  integration: 


_ _  C  3 _  /rp\ 

1  +  4™*At  r  Sr  'rD' 


4:1  JB  +  1  I 


1  +  4TTO*At 


(2.30) 


This  result  is  a  more  parallel  treatment  of  Ez  and  E  ,  treating  only  the  magnetic 
field  and  the  beam  current  Jg  as  source  terms.  The  exponential  treatment  has  the 
aovantage  of  eliminating  the  initial  value  ^  faster,  but  has  the  added  complica- 
tio„  of  ("*)(4i;  ft)  appearing,  which  can  lead  to  sign  changes  in  various  terms 
where  4:  ;t  is  not  small.  Note  that  in  Eq .  /.3  O'.  -B/  t  acts  directly  as  a 

source  of  E  ,  and  (rB)  acts  directly  as  a  source  of  E  ,  unlike  the  situation 
in  the  formulation  described  in  Section  2.3.2  above. 

The  attractiveness  of  this  approach  degrades  when  the  next  step  is 
examined.  Since  Eq .  (2.18)  cannot  be  written  in  terms  of  f  and  :  only,  it  is 
necessary  to  split  (2.30)  into  its  real  d  imaginary  parts.  In  principle  both 
Er  and  Ez  separately  from  (2.30)  can  be  inserted  in  Eq.  (2.24),  yielding  a  single 
equation  for  B.  However,  because  of  the  SB/St  term  in  the  expression  for  Ez  and 
the  spatial  derivative  in  the  expression  for  E  ,  the  resulting  equation  does  not 


& 


& 


involve  B  in  a  diagonal ly-domi nant  form,  and  thus  is  much  less  convenient  to  solve 
accurately.  In  addition  it  is  much  more  complicated  than  Eq.  (2.28)  because  of 
the  splitting  of  (1  +  4i:o*it)’^  into  real  and  imaginary  parts. 


A  variation  of  this  approach  was  actually  tried  prior  to  beginning  the 
present  study.  This  consisted  of  using  Er  from  Eq.  (2.30)  to  eliminate  Er  from 
Eq.  (2.24),  which  was  written  in  the  form 
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°u  Er 


(2.31) 


This  equation  was  solved  simultaneously  with  the  Ez  equation  from  (2.30)  for  E 
and  B.  Several  minor  variations  of  this  were  tried.  In  one  case  Eq.  (2.31)  was 
centered  as 


-zk+l 
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^°l  Er^k+% 


(2.32) 


Since  all  field  variables  were  to  be  defined  on  the  same  spatial  grid,  the  k+^ 
centering  of  the  Er  expression  from  Eq.  (2.30)  is  not  inconvenient  because  of 
the  appearance  of  the  spatial  derivative  of  B.  However,  the  Ez  equation  from 
(2.30)  is  awkward  to  center  with  Ez  and  B  on  the  same  grid.  The  most  direct  way 
would  involve  B  at  three  successive  grid  points.  This  was  avoided  by  writing 
tne  left  hand  side  of  (2.30)  as 


\ 

2 


+ 


Ezk+1 
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and  using  the  natural  centering  on  the  right  hand  side.  It  was  found  that  this 
method  of  solution  was  highly  unstable  and  totally  unreliable. 


Another  variation  tried  involved  writing  (2.31)  as 


<>  *  T" 


r)  -  E 


zk 


’  r .  h .  s .  '■ 


v  2. 33) 


(in  which  r.h.s.  is  a  function  of  B  only),  when  c  >  0,  and  the  analogous  thing 
with  centering  in  the  other  direction  when  c  <  0.  The  Ez  equation  from  (2.30) 
was  treated  as  described  above.  This  system  of  two  equations  gave  reasonable- 
looking  solutions  most  of  the  time,  provided  that  c  was  large,  although  it  some¬ 
times  went  unstable  even  in  that  situation.  For  cases  where  c  /c  «  1  it 
frequently  exhibited  unstable  behavior.  Thus  it  was  impossible  to  do  a  study  on 
the  transition  from  small  to  large  c  using  this  one  algorithm. 

Apart  from  the  occasionally  unstable  behavior,  this  algorithm  has  a 
more  serious  problem  associated  with  the  spatial  centering  of  Eq.  (2.33).  That 
spatial  centering  is  seriously  biased  if 

4 me  Pr  , 


It  was  not  always  possible  to  avoid  this  situation  with  a  reasonable  number  of 

smoothly-distributed  spatial  grid  points  extending  to  large  radius.  This  problem 

could  have  been  avoided  by  a  more  complicated  treatment  of  the  right  hand  side. 

The  basis  for  such  a  treatment  can  be  seen  by  writing  the  solution  to  Eq.  (2.31) 

in  the  form:  r 

A  k+i 

/  •'  r  \  //  \  (r-Pr)  \ 


Ezk+1  =  Ezk  exp(-  — )exp( 


in  which  r-dependcnce  of  has  been  ignored  for  simplicity  of  exposition.  Even 

the  simple  assumption  that  :  E  /j  varies  m  linearly  with  r  over  pr  would  allow 
a  formulation  in  which  the  obvious  bias  of  Eq.  (2.33)  is  avoided.  However,  the 
cost  involved  is  that  the  magnetic  field  at  as  many  as  four  points  rather  than 
two  would  appear  in  the  analogue  of  Eq.  (2.33)  because  of  the  unnatural  centering 
required  in  Eq.  (2.30)  for  E  . 


There  is  another  more  fundamental  objection  to  any  simple  implementation 
of  Eq.  (2.33).  For  sake  of  discussion  assume  that  it  is  further  developed  by 
assuming  o  E  /c  is  -  linear  over  Lr  .  Then  it  is  clear  that  when 


it  will  enforce  the  relationship 

(°iEA+l  "  KEA+1 

on  the  fields.  Thus  it  is  equivalent  to  forcing  0,  but  is  independent  of 
the  time-evolution  of  the  solution  and  is  therefore  not  physical.  The  resolution 
of  this  paradox  is  that  o|;  E  /o,  must  be  allowed  to  have  exponential  behavior  in  r. 
Only  then  can  a  physical ly-consi stent  treatment  be  developed  along  these  lines. 

The  bad  behavior  of  the  scheme  using  the  correctly  centered  Eq.  (2.32) 
can  be  identified  by  recalling  the  cancellation  that  occurred  in  going  from  Eq. 
(2.25)  to  (2.26).  If  the  fields  from  Eq.  (2.30)  were  both  inserted  in  the  right 
hand  side  of  Eq.  (2.24)  and  the  cancellation  done  analytically,  the  remaining 
equation  could  then  be  centered  as  in  Eq.  (2.32).  Although  this  procedure  has 
not  been  tried,  it  is  likely  that  much  of  the  bad  behavior  of  the  system  based 
on  Eq.  (2.32)  would  be  eliminated. 

As  in  the  method  discussed  in  Section  2.3.2,  there  remains  the  problem 
of  the  non-1 inearity  introduced  by  the  strong  dependence  of  o,  on  B.  All  of  the 
methods  discussed  in  the  present  section  2.3.3  involve  o  in  many  more  places 
and  in  more  complicated  ways  than  in  Eq.  (2.28).  No  attempt  has  been  made  to 
iterate  them  because  of  the  general  success  of  the  method  discussed  in  Section  2.3. 


2.4 


Computational  Results 


2.4.1  Introduction 

Two  different  SAI  codes,  HICHEM  and  HIGAP,  were  used  to  perform  the 
calculations  described  below.  Both  incorporated  the  same  basic  electromagnetic 
algorithm  which  was  described  in  Section  2.3.2.  The  tensor  cases  were  iterated 
until  the  magnetic  field  changed  by  less  than  1*  between  consecutive  iterations. 

The  simple  model  of  the  Hall  current  described  in  Section  2.2  was  used  for  all  the 
case  studies.  One  comparison  case  was  repeated  using  the  full  theory  as  described 
in  Ref.  1. 

The  HICHEM  code  (Ref.  4)  is  a  radially-resolved  air  chemistry  code 
specifically  designed  to  follow  the  air  chemistry  processes  generated  by  an  electron 
beam.  It  follows  the  non-equilibrium  time  evolution  of  up  to  60  chemical  species, 
and  can  account  for  radiative  energy  losses  when  these  are  important  to  the  overall 
energetics.  The  conductivity  calculation  is  also  described  in  Ref.  1;  it  incor¬ 
porates  a  very  detailed  model  for  electron  collisions  with  neutral  air  species, 
ions,  and  other  electrons,  and  computes  the  off-diagonal  (Hall)  component  of  the 
conductivity  tensor  in  a  consistent  way.  The  electromagnetic  field  equations  were 
solved  on  a  grid  of  290  radial  points  extending  to  500  Bennett  radii.  The  geometric 
grid  spacing  provided  18  points  inside  1  Bennett  radius,  and  58  points  inside  5 
Bennett  radii.  Except  as  noted  below,  all  the  HICHEM  calculations  used  the  same 
31  chemical  species  and  283  reactions.  The  cases  which  were  followed  out  to  50  ns 
used  37  species  and  337  reactions  in  order  to  account  for  substantial  radiative 
cooling  of  the  longer  pulses. 

The  HIGAP  code  (Ref.  5)  is  a  monopole  propagation  code  designed  to  study 
monopole  envelope  evolution  and  nose  erosion  for  axi-symmetric  beams.  The  con¬ 
ductivity  calculation  in  this  code  is  provided  by  the  BMCOND  (Ref.  6)  package  of 
subroutines.  This  package  uses  a  highly  simplified  model  of  the  chemistry  which 
has  been  calibrated  against  the  more  complicated  HICHEM  code.  Several  important 
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changes  in  these  calculations  were  required  to  handle  the  high  degree  of  ionization 
produced  by  high-current  beams  in  low  density  air.  The  electromagnetic  field 
equations  were  solved  on  a  radial  grid  of  250  points  extending  to  20  cm  from  the 
beam  axis. 

Since  the  HIGAP  code  (in  its  non-propagating  mode)  is  much  less  costly 
to  run  than  the  HICHEM  code,  it  was  used  to  investigate  some  more  unusual  or 
extreme  cases  for  sensitivity  to  Hall  effects. 

For  all  calculations  except  Case  6,  as  noted  below,  the  beam  profile  was 
assumed  to  be  Bennett.  The  time-dependence  of  the  beam  current  was  given  by 

I  =  IB  tanh(^)  .  (2.34) 

The  expanded  nose  of  the  beam  was  described  by  a  time-varying  Bennett  radius 
specified  by 


in  which  a  is  the  final  radius  as  t  -*•  «.  The  initial  radius  and  the  position 
o 

and  steepness  of  the  pinch-down  are  controlled  by  the  parameters  b,  t  ,  and  t^. 
2.4.2  Case  Studies 

Numerical  results  are  discussed  below  for  the  cases  described  in  Table  2- 
The  beam  current  in  each  case  was  100  kA.  These  cases  were  selected  to  provide  a 
wide  range  of  beam  parameters  and  channel  properties. 


The  basic  model  for  the  hot  channel  cases  had  a  central  temperature  of 
6000°K  with  the  density  chosen  to  give  pressure  equilibrium  at  1  atmosphere.  Three 
channel  widths  --  flat  out  to  1,  2,  or  3  Bennett  radii  of  the  pinched  beam  (as  in 
Eq .  2.35)--  were  investigated.  For  cases  3  -  6,  a  cold  uniform  medium  with  density 
(and  pressure)  =0.1  normal  was  used. 

Nose  parameters  were  chosen  to  simulate  a  variety  of  situations.  For 
Cases  la  and  2c  the  parameters  are  interpreted  as  providing  a  tailored  emittance 
which  decreases  significantly  along  the  beam.  Case  3b  represents  a  beam  with  a 
moderately  expanded  nose,  in  contrast  to  Case  3a  which  is  a  cylindrical  slug  of 
charge.  The  parameters  of  Case  4  represent  a  very  mature  nose  after  considerable 
erosion  of  the  beam-head  has  occurred. 

Case  5  was  chosen  to  investigate  the  effect  of  an  initial  distribution- 
of  electron  density  provided  by  an  annular  laser  beam.  This  was  expected  to  change 
the  conductivity  development  in  the  wings  and,  hence,  the  distribution  of  plasma 
current. 

Case  6  had  a  hollow  beam  current  distribution,  which  resulted  in 
significantly  different  conductivity  and  plasma  current  distribution  from  the 
conventional  Cases  1-4.  The  beam  current  density  on-axis  was  about  .075  of 
the  maximum,  which  occurred  at  about  0.5  Bennett  radii. 

In  addition  to  the  fully  self-consistent  calculations  from  the  HICHEM  and 
HIGAP  codes  (Cases  1  -6),  two  artifically  constrained  situations  were  studied  in 
order  to  assess  the  sensitivity  of  the  calculations  to  the  conductivity  profile. 
Conductivity  profiles  with  (3)  the  square  root  of  the  beam  profile  and  (t>)  the 
Bennett  profile  of  the  beam  were  imposed.  The  axial  conductivity  was  taken  to  be 
proportional  to  the  beam  current,  and  there  was  no  feedback  to  the  conductivity 
from  the  rest  of  the  calculation. 

All  the  cases  listed  in  Table  2-1  were  calculated  using  both  the  scalar 
and  tensor  versions  of  the  electromagnetic  algorithm.  Detailed  discussions  of 
the  individual  cases  are  given  in  Section  2.4.4  below,  preceded  bv  a  neneral  dis¬ 
cussion  of  results  in  Section  2.4.-. 


2.4.3 


Discussion  of  Results 


For  a  self-pinched  beam  the  force  (Hall  term)  exerts  an  outward 
force  on  the  backwards-flowing  J  plasma  current.  If  this  effect  reduces  the 
plasma  current  density  near  the  beam  axis,  the  net  current  density  and  local 
magnetic  field  increase.  This  in  turn  leads  to  an  increased  force  on  the 
remaining  plasma  current.  If  this  non-linear  response  is  effective  in  expelling 
a  significant  amount  of  the  plasma  current  from  the  vicinity  of  the  beam,  the 
beam  is  more  strongly  pinched  and  its  propagation  characteristics  are  improved. 

The  calculations  show  that  with  or  without  Hall  currents,  the  electro¬ 
magnetic  fields  rapidly  evolve  toward  the  situation  Jr  'v  0.  The  magnetic  force 
term  which  pushes  the  plasma  current  outward  is  then  balanced  by  the  radial 
electric  field.  Only  a  small  residual  J  remains,  generated  by  longer  time-scale 
evolution  of  the  magnetic  field  (due  to  conductivity  profile  evolution,  for 
example,  or  continued  evolution  of  the  beam  current  or  its  radial  distribution). 

In  none  of  the  self-consistent  cases  investigated  was  there  a  large 
scale  redistribution  of  the  plasma  current.  Although  ratios  o_/c|(  s  0.25  were 
achieved  late  in  Case  la,  for  example,  J  (r)  was  essentially  unchanged  from  the 
scalar  case.  This  result  is  consistent  with  the  simple  theory  discussed  in 
Section  2.2,  because  the  condition  0  -^0  was  reached  long  before  :  /:  reached  a 
large  value. 

The  constrained  model  with  conductivity  profile  specified  to  be  the 
square  root  of  the  Bennett  beam  profile  was  the  only  case  to  show  significant 
differences  between  the  scalar  and  tensor  cases.  The  magnetic  field  at  1  Bennett 
radius  increased  steadily  and  essentially  linearly  with  time  from  before  0.01  ns 
to  beyond  10  ns.  The  scalar  and  tensor  calculations  even  in  this  case  did  not 
differ  significantly  until  about  1  ns.  In  contrast,  the  constrained  model  with 
a  Bennett  conductivity  profile  had  a  net  current  which  was  a  factor  4  lower 
at  0.01  ns,  and  rose  at  a  much  slower  rate  initially.  By  10  ns  the  ratio 


o  /cr  reached  %  0.25,  but  the  difference  between  scalar  and  tensor  calculations 
was  negligible,  in  accordance  with  the  simple  theory.  For  all  the  self-consistent 
cases,  the  early  (before  1  ns)  behavior  of  •:  /cM  was  qualitatively  similar  to  that 
of  the  Bennett  profile  constrained  case;  it  increased  relatively  slowly  in  time, 
or  even  remained  roughly  constant. 

It  is  tempting  to  conclude  from  the  comparison  of  the  two  constrained 
cases  that  the  conductivity  profile  determines  whether  or  not  Hall  current  effects 
can  develop.  For  self-consistent  calculations,  several  processes  operate  to 
determine  the  conductivity  profile  evolution.  Initially,  in  a  cold  gas  the  ioniza¬ 
tion  is  produced  by  the  beam  and  has  a  beam  profile  --  Bennett  in  all  but  Case  6 
of  the  present  study.  As  the  Coulomb  Er  field  grows  rapidly,  it  breaks  down  the 
air  off-axis,  causing  a  temporary  (usually)  off-axis  peak  in  o.  Later  the  Ez 
field  rises  on-axis  and  dominates  the  air  breakdown.  The  peak  Ez  field  may  occur 
as  early  as  0.05  ns  if  a  broad  nose  is  not  present  on  the  pulse.  Usually  the  Ez 
field  brings  the  peak  of  the  conductivity  profile  back  to  the  axis  and  gives  it  a 
shape  roughly  Bennett  or  even  slightly  narrower.  As  the  conductivity  grows 
rapidly  near  the  axis  the  Er  field  begins  to  short  out  and  its  peak  moves  outward 
beyond  the  main  body  of  the  beam.  Conductivity  in  the  far  wings  rises  as  the  air 
is  broken  down  by  Er.  If  the  electron  production  on-axis  is  not  too  great,  N0 
begins  to  saturate  due  to  dissociative  recombination  on  molecular  ions.  In  this 
situation  a  profile  comparable  to  the  square  root  of  the  Bennett  can  develop  (but 
generally  at  t  -•  1  ns).  On  the  other  hand,  if  the  gas  on-axis  is  ionized  too 
quickly  for  recombi nation  to  keep  up  (or  molecular  ions  are  destroyed  by  the  high 
temperatures)  the  Spitzer  conductivity  regime  may  be  reached.  As  noted  in  Section 
2.2  above,  this  generally  sharpens  the  conductivity  profile  neartheaxis,  although 
broad  wings  will  persist.  This  entire  sequence  can  be  modified  by  providing  a 
broad  initial  hot  channel  (as  in  Cases  1  and  2),  an  expanded  nose,  a  long  taper  of 
the  beam  radius  due  to  a  variable  enittance,  or  by  a  non-Bennett  beam  profile 
(Case  6)  or  a  laser-prepared  channel  (Case  5).  A  very  wide  variety  of  profiles 
was  provided  by  the  self-consistent  cases  investigated,  but  in  no  instance  did 
the  qualitative  behavior  of  the  constrained  broad-profile  Case  7a  emerge.  The 


feedback  between  the  fields,  net  current,  and  conductivity  is  evidently  more 
important  than  any  of  the  modifications  of  conductivity  profiles  achieved  by 
choice  of  channel  parameters,  beam  parameters,  nose  parameters,  etc.,  in  the 
numerical  calculations.  For  high  current  beams  in  low-density  air  the  con¬ 
ductivity  evolves  rapidly  enough,  and  the  induced  Ez  is  strong  enough  to  allow 
very  nearly  complete  current  neutralization  to  be  achieved  early  in  the  pulse, 
in  the  spatial  region  occupied  by  the  beam  itself  rather  than  through  large 
return  currents  in  the  wings. 

As  noted  above,  the  constrained  Case  7a  did  not  develop  large  differences 
between  the  tensor  and  scalar  calculations  until  after  %  1  ns,  but  both  were 
qualitatively  different  in  magnetic  field  behavior  (as  measured  by  oJa{])  from 
Case  7b  and  the  self-consistent  cases.  This  suggests  that  it  may  be  possible 
to  identify,  from  scalar  calculations  only,  situations  in  which  significant 
tensor  effects  may  arise. 


2.4.4 


Details  of  Case  Studies 


Case  1:  2  mm  beam  into  a  low  density  hot  channel 

Several  different  channel  widths  were  investigated,  as  well  as  the 
effect  of  a  variable  emittance  along  the  beam.  Case  la  achieved  the  largest 
ratio  of  c  /o  .  The  radial  profiles  of  c,  and  o,  are  shown  in  Figure  2.1  at 
several  times  late  in  the  pulse,  where  o  has  grown  to  significant  values 
0.25  c,:  at  2  Bennett  radii  from  the  axis.  In  spite  of  the  large  o  /or,  the 
difference  in  pinch  function,  net  current,  Ez  filed,  etc.  between  tensor  and 
scalar  calculations  is  negligible.  Only  the  field  is  substantially  different 
from  its  value  in  the  scalar  case,  and  is  closely  given  by  the  condition  Jr^0. 

The  net  current  integrated  to  radius  r  (including  displacement  current) 
is  shown  in  Figure  2.2  as  a  function  of  radius  and  time.  The  specific  volume 
relative  to  sea  level  air  is  also  shown.  It  is  clear  that  the  plasma  current 

density  exceeds  the  beam  current  density  in  the  region  between  about  2 

Bennett  radii  and  the  channel  wall  at  4  Bennett  radii.  Even  as  late  as  50  ns 
into  the  pulse,  the  net  current  inside  4  Bennett  radii  is  only  2  kA,  although 

inside  1  Bennett  radius  the  net  current  is  about  16  kA  at  that  time.  In  the 

region  r  <  a,  the  beam  current  is  %  84%  neutralized  by  plasma  current. 

For  Case  lb  with  channel  walls  moved  inward  by  2  Bennett  radii,  tensor 

effects  are  weaker  still  because  the  plasma  current  flow  is  even  more  strongly 

confined  to  the  same  spatial  region  as  the  beam  current.  Although  an  expanded 

nose  and  emittance  tailoring  could  have  improved  the  situation,  it  is  unlikely 

they  would  have  resulted  in  significant  differences  between  the  tensor  and  scalar 

calculations  before  J  ->0. 
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Case  2:  5  mm  beam  into  three  different  hot  channels 

The  conductivity  is  shown  for  this  lower  beam  current  density  in 
Figure  2.3.  The  channel  profile  is  the  same  as  in  Case  la.  Figure  2.1.  In 
Case  2c  the  lower  current  density  results  in  less  strong  plasma  currents  near 
the  axis,  but  the  magnetic  field  in  the  current  distribution  is  weaker  in 
general  than  in  Case  la.  The  ratio  c  /a.,  barely  exceeds  .10  even  late  in  the 
pulse.  Again,  this  is  long  after  the  radial  plasma  current  has  approached  zero, 
so  there  is  essentially  no  difference  between  the  tensor  and  scalar  computational 
results  except  for  the  E  field.  As  expected,  the  narrow  channels  (Cases  2a  and 
2b)  were  less  favorable  to  the  development  of  a  strong  magnetic  field  because  of 
the  confinement  of  the  plasma  current  to  the  region  near  the  beam. 

Case  3:  2  mm  beam  into  uniform  density  =  0.1  of  sea  level  value 

In  some  ways  this  case  is  like  the  extreme  of  a  very  wide  channel,  but 
with  a  very  important  qualitative  difference.  Since  the  Case  1-like  channels  are 
all  chosen  to  be  in  pressure  equilibrium,  they  have  a  significant  initial 
conductivity  10  sec  in  their  low  density  (high  initial  temperature) 
regions.  In  Case  3,  however,  there  is  no  initial  conductivity  channel  because 
the  temperature  is  taken  to  be  uniform  at  288° < 

Conductivity  results  are  shown  in  Figure  2.4.  Because  the  beam  itself 
is  confined  to  a  narrow  region,  it  (and  the  accompanying  ohmic  heating)  generates 
its  own  conductivity  channel  which  is  initially  relatively  confined  within  a 
few  Bennet  radii.  Again,  the  plasma  currents  overlap  the  beam  current  sub¬ 
stantially  at  early  times,  preventing  it  from  establishing  a  strong  magnetic 
field  early  in  the  pulse.  Integrated  plasma  current  profiles  are  shown  as  a 
function  of  radius  at  several  times  during  the  pulse  in  Figure  2.5. 
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Figure  2.5.  Case  3a:  (2mm  Radius  No  Nose) 

Net  current  Integrated  to  Radius  r 
I  .  =  .005  r  B(r)  (kA) 


In  Case  3b  an  expanded  nose  was  added  to  the  pulse  of  Case  3a.  This 
resulted  in  somewhat  earlier  development  of  conductivity  in  the  region  outside 
the  beam,  but  was  still  insufficient  to  prevent  a  high-density  plasma  current 
from  cancelling  much  of  the  beam  current. 

Case  3a  was  repeated  using  the  full  tensor  calculation  of  c  rather  than 
the  simple  model.  In  general  this  gave  slightly  larger  values  of  o:,  but  still 
too  small  at  early  times  to  have  any  influence  on  the  electric  or  magnetic 
fields  before  -«•  0.  Out  to  ID  ns  the  non-cancellation  of  tensor  effects  in 
the  more  complex  theory  of  the  Hall  current  led  to  only  small  deviations  from  the 
earl ier  calculation. 


Case  4:  5  mm  into  0.1  normal  density  cold  air;  mature  nose 

This  case  models  a  very  mature  nose  such  as  might  exist  after  a  consider¬ 
able  amount  of  erosion  has  occurred.  The  full  100  kA  is  flowing  in  the  expanded 
region  before  the  pinch-point.  This  initial  broad  beam  generated  before  the  rapid 
pinch-down  an  electron  density  of  ^  5  x  10^  and  conductivity  v  5  x  10^  sec  1  - 
much  broader  than  the  final  beam  radius  of  5  mm.  It  was  a  further  attempt  to  get 
the  plasma  return  current  to  flow  outside  the  beam.  Some  redistribution  was 
achieved,  as  in  Case  3b.  However,  the  degree  of  current  neutralization  in  the 
core  of  the  pinched  beam  was  not  reduced  sufficiently  for  Hall  effects  to  be 
important  in  the  pinch-down  region.  To  some  extent  the  generation  of  initial 
conductivity  in  the  wings  is  also  counterproducti ve  to  the  generation  of  Hall 
effects  early  in  the  pulse.  This  same  conductivity  which  allows  the  plasma 
current  to  move  outside  also  shorts  out  the  Coulomb  field,  allowing  the  Jr  ■+•  0 
condition  to  be  achieved  more  quickly.  As  noted  before,  this  works  against  strong 
Hall  effects. 


Case  5:  5  mm  beam  in  0.1  normal  density  cold  air 


An  initial  electron  density  profile  given  by 
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was  used  to  simulate  an  annular  ionized  but  cold  channel  provided  by  a  laser. 

The  peak  N£  is  at  3  Bennett  radii  from  the  beam  axis,  and  drops  off  quickly 
inside  2  BR  and  outside  4  BR.  This  initial  seeding  with  electrons  provides  a 
conductivity  in  the  range  10-10  sec  ,  and  changes  the  breakdown  characteristics 
of  the  channel.  However,  it  did  not  result  in  plasma  current  redistribution 
sufficient  to  promote  significant  Hall  effects  early  in  the  pulse. 


The  electron  density  profile  at  0.05  ns  is  shown  in  Figure  2.6  for 
Case  5  (the  laser  channel)  and  Case  6  (the  hollow  beam),  along  with  a  comparison 
profile  for  the  same  type  of  case  with  no  channel  or  hollowing.  For  Case  5  the 
peak  electron  density  on-axis  was  already  three  orders  of  magnitude  greater  than 
that  in  the  initial  laser  channel  and  the  Er  field  had  already  reversed  in  sign 
out  to  ^  2  Bennett  radii.  The  peak  electron  density  and  conductivity  occurred 
at  %  0.3  Bennett  radii,  the  result  of  earlier  breakdown  by  the  radial  Coulomb 
field.  Essentially  no  net  current  flowed  in  the  region  between  0.5  and  1.2 
Bennett  radi i . 


Case  6:  5  mm  beam  into  0.1  normal  density  cold  air;  hollow  beam  profile 


The  radial  profile  of  beam  current  was  calculated  from 
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instead  of  the  usual  Bennett  formula.  Parameters  used  were  a  =  0.5,  b  =  0.25, 
=  0.32.  With  this  choice,  the  beam  current  had  a  hollowed  profile,  with  the 
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axial  current  density  less  than  1/10  the  peak  value  which  occurred  near  0.5 
Bennett  radii.  This  profile  function  was  chosen  only  as  a  crude  representation 
of  what  a  physical  hollow  beam  profile  might  be  like. 


The  distribution  of  plasma  current  was  found  to  remain  very  similar  to 
the  distribution  of  beam  current  up  to  about  1  ns,  with  a  high  degree  of  current 
neutralization.  By  1  ns  the  beam  current  density  had  reached  its  peak  value,  but 
the  plasma  current  density  on-axis  continued  to  grow  to  about  50%  higher  than  the 
beam  current  density;  the  magnetic  field  was  reversed  inside  %  0.25  Bennett  radii. 
The  conductivity  profile  remained  peaked  off-axis  through  the  10  ns  of  the  cal¬ 
culation.  The  net  current  integrated  out  to  large  radius  was  a  few  percent  higher 
than  that  for  Case  5  at  10  ns,  but  a  factor  of  about  2  lower  than  for  Case  4. 

As  in  those  cases,  there  was  no  significant  difference  between  the  tensor  and 
scalar  calculations. 

Case  7:  2  mm  into  uniform  low  density  channel;  artifical ly-imposed 

conductivity  profile 

The  purpose  of  this  calculation  was  to  assess  the  sensitivity  of  the 
electromagnetic  fields  to  an  artifical ly-imposed  conductivity  profile.  This  was 
done  by  decoupling  the  conductivity  calculation  completely  from  the  fields.  The 
beam  current  was  calculated  as  described  by  Eq.  (2.34)  above,  and  the  conductivity 
on-axis  was  simply  chosen  to  be  proportional  to  the  beam  current,  with  a  peak 
value  of  1014  sec'1  at  I  =  100  kA.  (a)  has  a  radial  profile  which  is  the 
square  root  of  the  Bennett  beam  profile  and  (b)  c  has  a  radial  profile  the  same 
as  the  beam  current.  The  ratio  -  /:  ,  and  /  itself,  were  computed  according  to 
the  simple  Hall  theory  described  above.  The  collision  frequency  was  chosen 
to  represent  a  uniform  low  density  region  with  .  v  0.036  of  normal  air  density. 

Tensor  and  scalar  results  for  the  two  different  profiles  are  shown  in 
Figure  2.7.  The  axial  electric  field  and  axial  plasma  current  density  are  shown. 
For  the  Bennett  profile  case,  the  beam  current  is  very  highly  neutralized;  the 
tensor  and  scalar  results  are  virtually  identical  out  to  20  ns.  On  the  other  hand, 
for  the  flatter  conductivity  profile,  a  significant  reduction  in  Ez  and  plasma 
current  occurred  in  the  tensor  case  beyond  1  ns. 
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2.5 


Cone! usions 


From  the  local  conductivity  studies  described  in  this  section,  several 
conclusions  emerge. 

(1)  In  order  for  Hall  effect  to  radially  redistribute  the  plsams 
currents  it  is  necessary  that  a  s  („ ,^/v  )  become  significant 
very  early  in  the  pulse  before  the  radial  Coulomb  field  reverses 
and  vanishes,  typically  within  the  first  few  tenths  of  a 
nanosecond.  If  this  condition  can  be  achieved,  on-axis  plasma 
current  decreases,  net  current  increases,  a  increases,  and  the 
effect  feeds  on  itself. 

(2)  At  low  gas  density,  v  decreases,  but  so  does  the  net  current- 
decreasing  ci,  Faster  avalanche  breakdown  or  ambient  ioniza¬ 
tion  result  in  Jr  tending  to  zero  earlier  in  the  pulse,  allow¬ 
ing  less  time  for  Hall  effect  to  act. 

(3)  The  strong  non-linearity  of  the  field  equations  -  including 
Hall  terms  -  requires  great  care  for  their  stable  solut  on. 

(4)  To  date,  we  have  numerically  studied  a  range  of  "realistic" 
cases,  with  conductivity  evaluated  sel f-consistently  with  the 
fields  but  always  assuming  local  conductivity  generation.  In 
all  these  cases,  the  race  is  lost;  current  neutralization  occurs 
so  fast  that  a  never  becomes  significant  until  well  into  the 
body  of  the  pulse  when  Hall  effects  essentially  cancel. 

Critical  to  all  these  results  is  the  conductivity  and  plasma  response  to  the 
beam  during  the  first  few  tenths  of  nanoseconds  -  a  time-scale  for  which  our 
local,  instantaneous  conductivity  modelling  is  not  appropriate,  especially  at 
low  densities.  More  adequate  non-local  conductivity  models  are  described  in 
the  next  section. 
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3.0 


NON-LOCAL  CONDUCTIVITY  EFFECTS 


3.1  Introduction 

This  section  summarizes  results  from  an  Initial,  simple  non-local, 
non-Ohmlc  air  chemistry  model  -  Including  Hall  effects  and  their  Impact  on 
beam  stabll Ity. 

In  Section  3.2  the  results  of  modelling  high  current  beams  In  low- 
density  air  are  discussed.  At  low-densities,  the  standard  simplifying  assump¬ 
tions  usually  employed  In  conductivity  modelling  no  longer  apply:  scalar  con¬ 
ductivity,  Ohm’s  law;  local-instantaneous  energy  deposition;  Maxwellian  dis¬ 
tributed  plasma  electrons;  no  delta  rays;  and  no  Inertial  effects.  These 
assumptions  are  not  made  In  the  present  model.  It  Is  concluded  that  Hall  cur¬ 
rents  do  play  a  significant  role  at  low  enough  densities  and  that  the  redis¬ 
tribution  of  plasma  current  can  result  In  a  significant  but  sudden  Increase  In 
the  magnetic  pinch  below  a  "critical”  air  density.  The  conductivity  model 
described  here  Is  a  slmplfled  version  of  the  LOCOND  model  described  elsewhere 
(Ref.  7)  and  has  been  developed  as  an  Intermediate  step  toward  a  truly  simple 
model  for  Inclusion  In  beam  stability  codes. 
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3.2 


High  Current  Beams  in  Low  Density  Air 


3.2.1  Introduction 

We  report  here  the  results  of  applying  a  phenomenological  low- 
density  chemistry  code  to  high  current  beams.  A  brief  account  of  the  model  has 
been  given  previously  (Refs.  8  and  9).  However,  substantial  modifications  were 
necessary  for  application  to  high  current  beams.  The  model  is  not  considered 
complete,  and  there  are  no  independent  calculations  in  the  same  parameter  regime 
to  which  it  can  be  compared;  thus.the  quantitative  results  are  tentative. 

Conclusions  based  on  the  calculations  are: 


1)  Below  a  certain  model -dependent  density  ^  .01  normal,  the 
electric  field  drives  a  bulk  runaway  which  changes  the  dis¬ 
tribution  of  plasma  currents.  The  net  result  is  a  sudden 
significant  increase  in  the  pinch  force  as  the  density  is 
decreased  below  the  critical  value. 

2)  If  Hall  currents  are  turned  off,  the  enhanced  pinch  may  be 
reduced  by  as  much  as  a  factor  of  two  at  one  Bennett  radius. 

3)  Significant  amounts  of  plasma  current  are  driven  by  the 
gradient  in  electron  pressure. 

4)  Significant  amounts  of  plasma  current  are  carried  by  the 
high-energy  non-Maxwell ian  part  of  the  electron  distribution. 

This  section  is  divided  into  three  main  topics:  (1)  a  discussion 
of  the  physics  requirements  for  low-density  calculations,  {l)  an  outline  of 
the  present  status  of  the  model,  and  (3)  a  discussion  of  the  computational 
results  for  high  current  beams. 


I 


3.2.2  Beam-Driven  Chemistry  in  the  Low-Density  Regime 

Simple  order-of-magni tude  argments  show  that  at  background  densities 
sufficiently  lower  than  normal  atmospheric  density,  several  key  assumptions 
built  into  standard  beam  chemistry  codes  are  not  justified.  Among  these  in¬ 
applicable  assumptions  are: 

1)  The  beam-initiated  cascade  results  in  local,  instantaneous 
production  of  electrons  and  ions. 

2)  Currents  associated  with  the  cascade  itself  can  be  neglected. 

3)  Plasma  currents  can  be  calculated  from  Ohm's  law. 

4)  The  electron  and  ion  densities  are  always  almost  equal,  so 
transport  effects  can  be  ignored  in  the  chemistry  calculations. 

5)  The  electron  distribution  remains  close  to  Maxwellian,  so 
that  the  high  energy  parts  of  the  distribution  are  no  more 
important  than  usual  in  determining  currents,  ionization 
rates,  etc. 

In  addition,  many  beam  chemistry  codes  use  electromagnetic  algorithms 
which  ignore  Hall  current  effects.  This  assumption  can  break  down  in  two  ways 
at  low  density:  (a)  the  momentum  transfer  frequency  goes  down  with  the  density, 
so  that  it  may  not  exceed  the  Larmor  frequency  by  a  wide  margin  as  in  full- 
density  air,  and  (b)  the  hi ghly-overpopulated  high  energy  tail  of  the  electron 
distribution  at  low  density  may  have  a  momentum  transfer  frequency  considerably 
lower  than  that  of  the  bulk  of  electrons,  and  thus  make  a  larger  contribution 
to  the  current. 

In  Section  2,  we  have  investigated  the  effects  of  Hall  currents  in 
full  and  reduced  density  air,  but  did  not  address  problem  (b)  above.  In  addi¬ 
tion,  we  have  made  assumptions  1  -  5  in  our  chemistry  codes.  We  have  recently 
developed  a  mul ti -energy-group  model  which  abandons  assumptions  1  -  5  and 
addresses  the  non-Maxwell ian  aspects  of  Hall  current  calculations  (Ref.  7). 

Recent  calculations  by  Yu  (Ref.  12)  have  confirmed  the  importance  of  Hall 
current  effects  for  ATA-like  beams  when  the  electron  energy  distribution  is 
treated  in  detail.  The  present  work  verifies  their  importance  for  hiqh  current 
beams . 
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3.2.3  Non-Local,  Non-Ohmic  Conductivity  Model 

The  ultimate  goal  of  the  development  program  is  to  produce  a  relatively 
simple  BMCOND-like  model  which  can  be  incorporated  into  propagation  codes.  This 
requirement  eliminates  the  possibility  of  doing  a  "first  principles"  calculation. 
The  approach  taken  is  to  develop  a  phenomenological  model  which  represents  the 
most  important  physical  processes  in  a  simple  way,  and  to  adjust  the  "free" 
parameters  associated  with  various  simplifying  assumptions  to  obtain  agreement 
with  more  fundamentally-based  calculations. 

The  model  described  below  considers  the  electron  distribution  broken 
into  three  energy  ranges.  The  lowest  energy  group  represents  basically  the  usual 
bulk  of  approximately-Maxwel 1 ian  electrons.  The  next  higher  energy  group 
represents  those  that  in  the  presence  of  an  electric  field  are  in  the  runaway 
regime  and  thus  behave  very  differently  from  the  bulk  of  the  electrons.  The 
third  group  represents  the  relativistic  particles  produced  directly  by  the  beam. 

The  organization  of  the  model  is  summarized  in  Fig.  3.1,  1n  which  the 
sources  and  sinks  of  particles  for  each  group  are  shown  schematically.  The 
three  energy  ranges  will  be  referred  to  by  the  names  low-group,  high-group, 
and  i  -group, in  order  of  increasing  energy.  The  beam  particles  themselves 
constitute  a  fourth  group  which  is  treated  in  the  usual  way  (e.g.,  no  energy 
straggl ing ) . 

The  beam  particles  collisionally  produce  secondaries  directly  in  each 
of  the  three  groups,  according  to  a  Moller  distribution  (modified  at  low  energy). 
The  maximum  5-ray  energy  is  one-half  the  beam  energy,  since  the  higher  energy 
particle  emerging  from  a  primary  interaction  remains  associated  with  the  beam. 
Collisions  by  beam  particles  are  the  only  source  of  electrons  in  the  5 -group. 

A  more  detailed  discussion  of  the  5-ray  model  is  given  below  in  Section  3.2. 3.1. 
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The  high-energy  group  extends  downward  from  1  MeV  to  a  variable-energy  boundary 
which  at  present  goes  no  higher  than  10  KeV.  It  is  fed  mainly  by  primary  beam 
collisions  and  by  the  collisional  degradation  of  the  5-group  through  its  1  MeV 
lower  boundary.  Under  some  conditions,  electrons  are  also  transferred  up  from 
the  low  group. 

Both  the  high  and  low  energy  electron  groups  are  represented  by  fluid 
equations  for  particle  density,  momentum  and  energy.  The  equations  of  motion 
provide  the  means  for  abandoning  the  ohmic  representation  of  plasma  current 
flow.  The  fluid  equations  are  solved  on  an  Tulerian  grid  which  is  the  same 
grid  used  for  the  electromagnetic  field  calculation.  It  is  assumed  that  every¬ 
thing  is  a  function  of  the  "retarded"  time  variable  Z  =  t  -  z/c  only,  so  that 
the  only  independent  variables  are  r,  the  radial  position,  and  £,  the  distance 
back  from  the  pulse  head  in  seconds.  The  electromagnetic  fields  are  obtained 
using  the  algorithm  described  in  Section  2. 

In  principle  it  is  necessary  to  solve  the  fluid  equations  and  5 -group 
equations  simultaneously  with  the  electromagnetic  field  equations.  Since  the 
Maxwell  equations  themselves  are  explicitly  non-linear  when  the  Hall  current 
terms  are  included  and  the  fluid  equations  are  also  non-linear,  an  iterative 
procedure  is  required.  In  order  to  maintain  maximum  modularity  under  these 
ci rcumstances ,  the  procedure  adopted  involved  solving  first  a  linearized  set 
of  fluid  equations,  followed  by  the  field  equations,  and  then  iterating  to 
convergence. 

In  addition,  the  algorithm  has  been  designed  to  allow  differential 
comparisons  between  the  three-group  model  and  a  model  with  the  high  energy 
and  5-groups  eliminated  and  with  the  standard  local-instantaneous  approximations 
made.  However,  this  latter  "one-group  model"  still  differs  from  standard 
chemistry  codes  in  that  it  solves  an  equation  of  motion  for  the  one  group 
rather  than  using  Ohm's  law.  It  is  also  possible  to  turn  off  the  pressure 
terms  in  the  equations  of  motion  and  the  Hall  current  effects 
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3.2.3. 1  Delta-Ray  Group 


4. 


In  the  present  discussion,  5-rays  will  refer  to  the  part  of  the  beam- 
initiated  cascade  which  has  energy  greater  than  1  MeV.  These  particles  are 
represented  py  a  model  essentially  identical  to  that  developed  by  Johnston 
(Refs.  7,  11).  As  applied  here,  the  model  has  two  main  deficiencies:  (1)  it 
does  not  provide  detailed  information  on  the  radial  distribution  and  (2)  it 
does  not  include  the  detailed  effects  of  electric  or  magnetic  fields  on  the 
-  particles.  With  the  assumption  that  everything  depends  on  £  e  t  -  z/c, 
the  model  gives  an  integral  expression  for  the  total  6-ray  current  as  a  function 
of  distance  •;  back  from  the  beam  head.  (Particle  energy  distribution  information 
is  available  but  not  used  at  present.) 

It  is  assumed  that  most  of  the  particles  produced  as  secondaries  by 
the  5-rays  have  very  low  energy.  This  is  consistent  with  the  continuous  slowing 
down  approximation  which  forms  the  basis  of  the  model.  In  the  present  model, 
as  the  5-rays  lose  energy,  they  produce  low  energy  electrons  at  a  rate  of  one 
particle  per  33.73  eV  of  energy  lost;  these  particles  immediately  enter  the  low- 
energy  plasma  electron  group  of  the  three-group  model.  The  5-rays  themselves 
degrade  in  energy  due  to  collisions.  They  reach  the  cutoff  of  1  MeV  at  a  rate 
given  by  the  Johnston  model,  and  at  this  point  they  are  added  to  the  high  energy 
electron  group. 

The  radial  profile  as  a  function  of  distance  back  from  the  beam  head 
is  estimated  in  a  crude  way,  based  on  two  pieces  of  information:  (1)  the  i's 
are  produced  with  the  beam  profile  and  (2)  they  evolve  from  that  profile  to 
form  a  halo  with  radial  dimension  estimated  by  Johnston  (Ref.  11)  to  be  -  4 
Bennett  radii  about  a  Bennett-prof i 1 e  beam  (for  particles  above  1  MeV).  It 
is  assumed  that  at  a  given  the  profile  towards  which  the  I's  evolve  is  a 
weighted  average  of  the  beam  profile  and  a  specified  halo  profile.  The  relative 
weights  are  determined  from  the  rate  of  local  production  as  compared  to  the 
total  instantaneous  number  density,  and  the  rate  of  evolution  towards  tne  nalo 
profile.  Tnis  a'su.tion  is  used  in  tne  following  interpolation  formula  for 
tne  radial  profile  f ,  r,  }  in  terms  of  f  r ,  : : 


f(r,  C+AC)  =  (f(r.c)  +  ACF)/(1  +  vdA£). 

The  function  F  is  given  by 

F  =  h  +  g  v), 

where  S  is  the  local  rate  of  5-ray  production  by  the  beam,  N  is  the  total  number 
of  5-rays  present,  and  h  is  the  beam  profile,  From  which  the  5-rays  are  produced. 
The  profile  g  is  the  assumed  halo  profile,  and  v  is  a  characteristic  rate  at 
which  the  5 ' s  evolve  towards  the  profile  g.  The  rate  vd  is  set  equal  to  S/N  +  v. 
With  this  choice,  f(r,£+A£)  is  normalized  to  unity  if  h,  g,  and  f(r,£)  are. 

The  rate  v  is  taken  to  be  an  average  radial  velocity  of  the  5 ' s  divided  by  the 
assumed  halo  characteristic  radius. 

3. 2. 3. 2  The  High-Energy  Group 

The  schematic  forms  of  the  fluid  equations  for  the  high-group  are: 

+  V*(Nv)  +  VjN  =  +  ^5  +  S|_  (3.1) 

4  +  v.(37)  ♦  I  VP  ♦  V  -  L  *  K.  ♦  K,  ♦  ^  (E  +  f  *B)  (3.2) 

dt  m 

(j  Nmv2  +  j  P)  +  V*[v(^  Nmv2  +  P)]+  7*(vP) 

+  Vj(  \  Nmv2  +  |  P)  =  L  +  NeE*v  (3.3) 

The  terms  SD,  KD,  etc.  are  the  sources  due  to  the  beam,  the  5-group,  and  upward 

D  D 

transfer  from  the  low-group,  respectively.  It  is  assumed  that  all  collisional 
ionizations  by  high-group  electrons  produce  a  low  energy  secondary  which  goes 
directly  into  the  low  group.  The  energy  equation  (3.3)  has  been  written  with 
the  total  energy  split  into  a  drift  part  and  a  thermal  part,  represented  by  a 
pressure  P.  The  term  L  includes  energy  deposited  by  the  beam,  energy  brought  in 
by  transfer  of  particles  from  the  low-group  and  i>-group,  and  energy  losses  due 
to  i  on  i  rati  or..  •  n*-r  :/  loss  due  to  loss  of  particles  from  the  high  group  is 
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represented  by  the  Vj  term  on  the  left  side  of  the  equation.  It  has  been  assumed 
that  the  particles  lost  take  with  them  the  average  energy  of  the  high  group.  The 
momentum  loss  rate  includes  the  effect  of  particle  losses,  momentum  loss  due 
to  ionizing  collisions,  and  momentum  loss  due  to  elastic  collisions. 

The  fluid  equations  are  solved  as  a  finite  difference  system  on  the  Eulerian 
electromagnetic  field  grid,.  The  outer  boundary  condition  allows  an  outflow  of 
particles.  No  explicit  assumption  is  made  concerning  the  distribution  function 
of  the  high  group  electrons  except  at  £  =  0+,  when  it  is  dominated  by  the  Moller 
distribution  from  primary  beam  particle  interactions.  The  average  total  energy 
and  z-momentum  are  then  used  to  give  an  initial  "temperature"  T,  defined  by  P  =  NT, 
where  P  is  the  pressure  from  Eq.  (3.3).  At  all  subsequent  times,  the  temperature, 
particle  density,  and  drift  velocity  are  obtained  from  the  fluid  equations. 

For  simplicity,  ionization  and  momentum  transfer  rates  are  evaluated 
at  the  mean  particle  energy.  An  analytic  formula  given  by  Briggs  and  Yu  (Ref.  12) 
is  used  for  the  ionization  cross  section.  The  momentum  transfer  frequency  vm  is 
given  by  the  approximate  formula 

-  =  (1.08  x  10‘5  N  £'2)/(178.89  +  £3/2) 

m  g 

+  (1.46  x  10"6  N+  log  A )/Tgq2-j y  (3*4) 

where ^ is  tne  total  electron  energy  in  eV.  Tne  first  term  represents  the 

effect  of  collisions  with  neutral  particles  and  the  second  represents  Coulomb 

2 

collisions.  The  "equivalent"  temperature  is  defined  as  Tg^u^v  =  T  +  1/3  mv  , 
log  /.  is  the  usual  Coulomb  logarithm,  Ng  is  the  total  density  of  neutral  atoms 
and  molecules,  and  N+  is  the  total  density  of  positive  ions. 

The  boundary  between  the  high  and  low  energy  groups  is  set  by  finding 
tne  solution  of  the  equation 

e  E  =  m  v  (35) 

m 
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corresponding  to  the  high  energy  side  of  the  peak  in  v  .  The  average  velocity 

m  ip  ^  2 

v  is  related  to  the  total  energy &  appearing  in  Eq.  (2.4)  by£  =  h  m  v  . 
However,  the  boundary  energy  is  not  allowed  to  exceed  10  KeV  even  when  the 
total  electric  field  E  is  very  small. 


Another  special  case  occurs  when  the  field  is  so  strong  that  there  is 
no  physical  solution  to  (3.5);  this  corresponds  to  a  bulk  runaway,  in  which  the 
field  is  capable  of  accelerating  the  electrons  through  the  peak  of  the  momentum 
transfer  cross  section.  In  this  case,  the  choice  of  the  lower  boundary  of  the 
high  group  is  somewhat  arbitrary.  The  lower  limit  to  be  used  for  the  boundary 
during  a  strong  runaway  condition  can  be  specified  as  input  data.  Usually  it 
is  taken  to  be  a  few  volts,  so  that  essentially  all  the  electrons  produced 
directly  by  the  beam  go  into  the  low  end  of  the  high  energy  group.  The  transfer 
of  electrons  from  the  low  to  the  high  group  is  discussed  below. 


3.2.3. 3  The  Low-Energy  Group 

The  equations  for  the  low  energy  group  are  similar  in  form  to  those 
given  above  for  the  high  group  with  the  exception  that  they  contain  loss  terms 
due  to  recombination  and  to  transfers  to  the  high  energy  group,  and  input  terms 
due  to  transfers  from  the  high  group  and  to  collisional  ionization  by  high 

group  electrons.  The  main  difference  in  the  treatment  of  the  low  group  is  in 

2 

tne  chemistry  detail.  At  present,  the  abundance  of  N^,  C^,  N,  N(  D),  0,  and  a 

composite  representati ve  of  the  triplet  states  of  are  followed  explicitly  by 
differential  equations.  The  system  of  equations  used  is  essentially  the  same 
as  in  the  chemistry  code  BMCOND  (Ref.  13).  This  degree  of  complexity  was  found 
necessary  to  obtain  reasonable  agreement  with  the  comprehensive  code  HICHEM  at 
low  air  densities  for  high  current  cal cul ations . 
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Ionization  rates  are  given  by  tables  identical  to  those  in  the  HICHEM 
code  up  to  500  volts  "equivalent"  temperature  (TeqU^v  =  Tg  +  1/3  m  v  ).  The 
collisional  excitation  rates  from  the  BMCOND  model  are  used.  Ionization  and 
excitation  rates  due  to  the  beam  are  computed  as  in  the  HICHEM  code,  but  with 
appropriate  adjustment  for  the  differences  between  time-delayed  deposition  model 
used  here  and  the  instantaneous  deposition  model  in  HICHEM  and  BMCOND. 

It  is  very  important  to  distinguish  between  the  positive  ion  density 
and  the  electron  density,  since  electron  transport  effects  can  be  very  large; 
thus  a  separate  differential  equation  for  the  total  rate  of  positive  ion  produc¬ 
tion  is  integrated.  In  addition  it  was  found  necessary  to  calculate  a  vibrational 
temperature  because  of  significant  sensitivity  of  computational  results  to  the 
dissociative  recombination  rate;  the  HICHEM  treatment  was  adopted  for  this 
calculation.  The  momentum  transfer  frequency  used  in  BMCOND  is  used  for  low 
energies,  with  a  smooth  transition  to  the  analytic  formula  (Eq.  3.4)  used  for 
the  high  group.  The  Coulomb  term  of  Eq.  (3.4)  is  used  throughout. 

Reassignment  of  electrons  from  the  low  group  to  the  high  group  is  a 
relatively  arbitrary  procedure.  The  goal  is  to  remove  electrons  from  a  presumed 
high-energy  non-Maxwel 1 ian  tail  of  the  low  group  and  put  them  into  the  high 
group,  which  hopefully  represents  better  their  contribution  to  currents,  ioniza¬ 
tion  rates,  etc.  than  does  the  assumed  Maxwellian  bulk  of  the  low  group  electrons. 
However,  there  is  no  simple  model  of  the  super-thermal  tail  as  generated  by 
strong  electric  fields  which  vary  rapidly  in  both  time  and  space.  Several 
different  ad  hoc  procedures  have  been  tried,  but  none  is  especially  defensible 
in  detail . 

Qualitatively,  what  happens  in  a  weak  field  (or  high  ambient  density 
gas)  is  not  very  sensitive  to  the  details  of  the  transfer  rate.  Under  strong 
runaway  conditions,  the  entire  low  group  is  accelerated  to  high  drift  energy 
and  heated  to  temperatures  -  kilovolts  on  a  very  short  timescale.  During  this 
time,  the  distinction  between  the  two  groups  is  not  very  meaningful,  and  again 
the  behavior  of  the  bulk  plasma  is  not  too  sensitive  to  the  details  of  the  transfet 
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rate.  However,  after  the  fields  which  precipitate  a  strong  runaway  die  out, 
the  distribution  of  electrons  between  the  two  groups  can  be  important;  this 
does  depend  strongly  on  the  transfer  rate.  Similarly,  if  conditions  for  a 
strong  runaway  are  just  barely  achieved  for  only  a  short  time,  there  may  be 
some  sensitivity  to  the  transfer  rate.  These  problems  are  continually  being 
studied,  both  in  the  context  of  the  model  described  here  and  through  comparison 
with  other,  more  comprehensive  models  being  developed. 


The  transfer  rate  presently  used  is  calculated  as  follows.  A  velocity 

v  isdefined  by  v  =  +  v  +  .5  V, ,  where  the  +  sign  is  taken  if  the  radial  drift 
o  o  -  r  l 

velocity  is  parallel  to  Er,  and  a  negative  sign  is  taken  if  it  is  opposed.  The 
thermal  velocity  term  is  added  to  account  for  the  fact  that, at  high  enough 
temperature,  a  substantial  number  of  electrons  may  be  able  to  run  away  even 
if  the  bulk  drift  is  not  large  (or  even  parallel  to  the  Er  field).  The  transfer 
rate  of  electrons  is  then  given  by 


v  -v 

SL  =  min( .5,  .lNh/Ng)  min(l.,  exp(-~y--c) )/AC , 


(3.6) 


in  which  vc  is  the  velocity  corresponding  to  the  energy  boundary  between  the 
groups,  Nh  is  the  local  density  of  high-group  electrons,  Ng  is  the  local  density 
of  low-group  electrons,  and  ££  is  the  proposed  timestep.  The  purpose  of  the  first 
minimum  function  is  to  assure  that  the  high-group  density  is  not  changed  by  more 
tnan  10  during  the  step.  In  actual  computations,  the  transfer  rate  is  usually 
very  low,  or  else  as  high  as  permitted  by  the  first  factor  of  (3.6),  even  for 
exceedingly  small  timesteps. 


The  momentum  transferred  with  the  particle  is  taken  to  be  parallel 
to  its  total  drift  velocity,  with  magnitude  given  by  the  larger  of  the  velocities 
v  and  v  .  The  energy  transferred  is  the  average  energy  per  particle  of  the 
low  group,  plus  an  additional  amount  arranged  to  come  (by  suitable  terms  in 
the  equations  of  motion)  entirely  from  the  drift  energy  of  the  low  group,  to 
make  up  the  tot.-'l  drift  energy  of  the  electron  injected  into  the  high  group. 
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The  net  result  is  to  decrease  the  drift  energy  of  the  low  group  but  leave  its 
temperature  unchanged.  The  effect  of  the  energy  and  momentum  transfers  on  the 
high  group  depends  on  its  energy  and  momentum  at  the  time  of  the  transfer. 

Other  procedures  for  transferring  electrons  to  the  high  group  are  being  evaluated. 

3.2.4  Calibration  in  the  High  Current  Regime 

At  present  there  are  no  comprehensive,  first-principle  calculations 
to  provide  detailed  guidance  in  the  development  of  the  three-group  model  for 
high  current  calculations.  Yu  (Ref.  10)  has  presented  Bol tzmann-code  calcula¬ 
tions  for  ATA-like  parameters,  but  even  these  calculations  incorporate  assump¬ 
tions  which  are  not  always  justified  in  the  context  of  high  power  beam  plasmas. 

The  purpose  of  the  calibration  runs  discussed  below  is  to  check  that  the 
relatively  simple  model  described  in  Section  3.2.3  agrees  reasonably  well  with 
detailed  chemistry  codes  in  regimes  where  the  assumptions  of  those  codes  (see 
Section  3.2.2)  are  not  thought  to  be  seriously  in  error.  Results  of  calibration 
at  10  kA  have  been  presented  elsewhere  (Ref.  0).  A  similar  comparison  with  the 
HICHEM  (Ref.  4)  code  at  100  kA  and  gas  density  0.1  normal  is  given  below.  The 
rise  time  is  5  ns  and  the  Bennett  radius  is  0.5  cm. 

The  comparison  of  electron  densities  on-axis  and  at  one  Bennett  radius 

is  shown  in  Figure  3.2.  The  agreement  is  very  satisfactory.  The  effect  on  the 

electron  density  due  to  the  time  delay  for  the  three-group  model  is  not  large 

after  2  or  3  ns.  The  electron  density  on  axis  is  mainly  determined  by  the 

close  balance  between  collisional  ionization  and  dissociative  recombination 

after  a  few  ns.  The  recombination  rate  itself  decreases  significantly  as  the 

molecules  (and  molecular  ions)  are  depleted,  and  as  the  vibrational  temperature 

increases.  After  about  6  ns,  the  continued  increase  in  N  is  determined  mainly 

e 

by  the  decrease  in  the  recombination  rate,  with  a  large  part  attributable  to 
the  vibrational  temperature  dependence. 


A  similar  comparison  of  electron  temperatures  and  total  currents  is 
shown  in  Figure  3.3.  The  temperature  agreement  is  quite  good  from  1  ns  on. 

At  0.1  ns,  HICHEM  gives  Tg  %  18  volts,  compared  to  the  12  or  13  volts  given 
by  the  simple  models.  At  such  early  times  the  ohmic  approximation  used  in 
HICHEM  is  not  very  good,  so  the  ohmic  heating  rate  is  likely  to  be  incorrect; 
thus  it  is  not  clear  which  values  are  closer  to  the  truth. 

As  in  the  10  kA  comparison  (Ref.  3),  the  agreement  between  HICHEM 
and  the  simple  models  is  not  as  good  for  the  net  current  and  effective  current. 
This  difference  is  mainly  due  to  the  difficulty  in  matching  the  momentum 
transfer  frequency  calculated  in  HICHEM  by  a  very  simple  formula.  However, 
the  agreement  in  currents  is  still  acceptable,  particularly  since  the  most 
interesting  sensitivities  described  below  develop  before  1  nsec. 

Because  the  results  of  the  HICHEM  code  are  suspect  below  p/pQ  ^  0.1, 
especially  for  ?  <  few  ns,  it  cannot  be  used  for  calibration  comparisons  at 
lower  densities.  However,  many  of  the  results  presented  below  are  differential 
comparisons  or  sensitivity  studies,  and  thus  can  provide  useful  information  in 
spite  of  uncertainties  in  the  quantitative  results. 

3.2.5  Computational  Results 

3. 2. 5.1  Introduction 

The  beam  parameters  for  the  calculations  described  below  are: 

Current  =  100  kA 

Rise  time  =  5  ns 

Bennett  radius  =  0.5  cm 

Energy  =  10  MeV 
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The  discussion  below  is  divided  into  two  main  sections.  The  sen¬ 
sitivities  to  gas  density,  Hall  currents,  and  pressure  terms  are  described 
briefly  in  Section  3.2. 5.2,  and  the  remaining  sections  discuss  the  calculationa 
in  more  detail. 

The  main  conclusions  have  been  stated  earlier  (Section  3.2.1).  In 
brief,  it  has  been  found  that  the  pinch  force  is  strongly  dependent  on  the  gas 
density  below  some  threshold  value,  and  significantly  dependent  on  the  presence 
of  Hall  currents  and  pressure  terms.  The  plasma  current  and  electron  density 
radial  distributions  change  dramatically  over  a  very  small  range  of  density 
ratio  near  p/pQ  =  .01,  resulting  in  an  increase  of  pinch  force  by  a  factor 
of  3  or  more  for  a  density  change  of  only  n,  20%. 

3. 2. 5. 2  Sensitivity  Results 

Density  Sensitivity  of  Pinch  Force 

The  dependence  of  Ie^(r)  on  p/pQ  is  shown  in  Figure  3.4  for  the 
three-group  model,  and  in  Figure  3.5  for  the  one-group  model.  The  sharp  onset 
of  current  enhancement  begins  in  the  range  p/pq  n,  .0085  -  .01  for  the  three- 
group  model,  and  between  .005  -  .008  for  the  one-group  model.  In  both  cases, 
the  current  enhancement  is  a  factor  3.5  or  greater,  compared  with  a  base 
level  at  ./.  01.  These  large  effects  set  in  when  the  plasma  electrons  can 

be  sustained  in  a  state  of  bulk  runaway  for  several  tenths  of  a  nanosecond. 


The  sharp  density  threshold  results  from  the  sensitivity  of  the  charge 
neutralization  process  to  gas  density.  Several  effects  contribute.  The  beam 
production  of  positive  ions  is  directly  proportional  to  the  density;  the 
electrons  produced  by  the  beam  provide  "see  for  the  avalanche  process  of 
plasma  electron  collisions  The  peak  e-folding  rate  of  the  avalanche  is  also 
proportional  to  the  density.  These  two  processes  clearly  delay  the  charge 
neutralization  as  the  gas  density  is  decreased  and  allow  higher  radial  electric 
fields  to  develop  (assuming  the  rise  time  is  not  ■■  1  ns).  In  order  for 


neutralization  to  occur,  plasma  electrons  produced  by  various  processes  must 
move  out  of  the  spatial  region  occupied  by  the  beam,  leaving  the  positive  ions 
to  cancel  the  charge  of  the  bear  particles.  Higher  fields  help  move  the 
electrons  out  quickly,  but  that  process  dilutes  very  considerably  the  avalanche 
ionization  and  thus  slo.-.:  f  t  production  of  the  needed  positive  ions  close  to 
the  bear  a*is.  In  tne  trj-  IT-  regime,  avalanche  is  relatively  unimportant,  and 
trie  accumulation  ratt  c  •  L-.  .s  -  generated  positive  ions  determines  the  neutral¬ 
ization  time  and  the  pea*  fields. 

The  difference  in  density  threshold  between  the  two  models  is  caused 
by  (1)  the  time-delay  in  ionization  in  the  three-group  model  which  reduces  the 
effective  bean;  ionization  rate  by  more  than  a  factor  of  two;  and  (2)  the  effects 
of  the  non-Maxwel 1 ian  high-energy  group  on  the  ionization  rate  and  on  the 
movement  of  plasma  electrons  away  from  the  beam.  It  is  difficult  to  assess 
these  separately  because  beam  production  of  the  high-energy  and  5 -ray  groups 
(which  degrade  relatively  slowly  at  low  densities)  is  the  cause  of  the  time  ■ 
delay.  Since  qual i tati vely-simi lar  results  occur  whether  or  not  the  high- 
energy  group  is  included,  it  seems  that  the  time-delay  is  probably  most  directly 
responsible. 

Sensitivity  to  Hall  Currents  and  Electron  Pressure 

The  radial  pressure  gradient  in  the  equations  of  motion  for  the  high 
and  low  energy  electrons  can  drive  currents  both  radially  and  in  the  z-direction 
due  to  the  magnetic  part  of  the  Lorentz  force  (Hall  effect).  This  works  in 
conjunction  with  the  electrically  driven  currents.  The  result  of  deleting 
either  the  pressure  terms  or  all  Hall  effects  (by  zeroing  the  magnetic  force 
on  the  plasma  electrons)  is  shown  in  Figure  3.6.  Clearly,  the  Hall  currents 
have  the  largest  effect,  but  the  pressure  terms  are  not  negligible  either. 

Note  that  these  comparisons  show  the  cumulative  effect  of  removing  the  terms  - 
for  the  entire  calculation,  and  not  simply  the  contributions  to  the  current 
at  the  time  shown.  A  similar  calculation  was  done  with  only  the  high  group 
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pressure  turned  off.  The  effect  at  £  =  1  ns  was  considerably  smaller  than  shown 
in  the  figure.  Low-group  pressure  contributes  dominantly  during  the  bulk 
runaway  because  the  electrons  reach  temperatures  in  excess  of  a  kilovolt  and 
their  number  density  is  large.  However,  it  will  be  shown  in  later  discussion 
that  the  high-group  pressure-driven  Hall  current  is  important. 

3. 2. 5. 3  Comparison  of  1-Group  and  3-Group  Calculations  at  p/p  =  .008 

The  purpose  of  the  following  discussion  is  to  provide  a  more  detailed 
description  of  the  phenomena  which  lead  to  the  enhanced  pinch  force.  At 
p/po  =  .008,  the  3-group  model  shows  a  very  strong  effect,  whereas  the  1-group 
model  shows  very  little  because  its  density  threshold  is  somewhat  lower. 

Effects  on  Plasma  Current  Distribution 

The  net  current  (including  displacement)  integrated  out  to  radius  r 
is  given  by  Inet  =  .005  r  B(r),  where  I  is  in  kA,  B  is  in  gauss,  and  r  is  in  cm. 
The  effective  current,  which  measures  the  pinch  force,  is  the  beam-profile- 
weighted  average  of  .01  r  (B(r)  -  E  (r)).  The  value  of  the  effective  current 
integral  taken  to  radius  r,  and  the  net  current,  are  shown  in  Figures  3.7(a)  -  (h) 
at  various  distances  from  the  pulse  head.  Large  differences  are  apparent  by 
1  ns  and  persist  to  10  ns. 

At  0.1  ns,  the  net  currents  for  the  two  calculations  are  similar, 
but  the  effective  current  is  substantially  weaker  for  the  3-group  case  because 
the  radial  electric  field  is  higher  (due  to  slower  charge  neutralization). 

At  1  ns  the  effective  current  in  the  3-group  model  is  n,  28%  of  the 
beam  current  at  that  time,  whereas  in  the  1-group  model  it  is  only  'v  7.5%. 

The  net  current  profiles  imply  a  very  much  broader  plasma  current  distribution 
for  the  3-group  model,  although  net  currents  inside  5  cm  radius  (10  Bennett 
radii)  are  ,  .5  to  1  kA  in  both  cases  even  at  10  ns  into  the  pulse  where  the 
beam  current  is  almost  100  kA. 
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A  comparison  of  the  plasma  current  (Jz)  radial  profiles  at  1  ns  is 
shown  in  Figure  3.8.  The  breakdown  of  the  3-group  profile  into  contributions 
from  the  high  and  low-energy  electron  groups  is  discussed  below. 

Electron  and  Ion  Density  Profiles 

Major  qualitative  differences  in  the  evolution  of  the  electron  density 
are  apparent  in  Figures  3.9(a)  -  (d).  The  1-group  calculation  shown  in 
Figure  3.9(a)  has  the  usual  peak  slightly  off  axis  due  to  the  Er~initiated 
avalanche  and  is  not  unusual.  The  3-group  calculation  for  the  low  group 
(Figure  3.9(b))  is  much  more  interesting.  It  shows  a  large  peak  well  off 
axis  and  much  higher  density  at  large  radius.  The  high-group  density 
(Figure  3.9(c))  has  features  which  invite  interpretation  as  propagating  sound 
waves.  The  ridge  which  appears  earliest  in  time  is  associated  with  the  rise 
of  the  Er  field  and  the  second  begins  on  the  axis  near  the  peak  of  Ez-  There 

is  a  hint  of  the  first  ridge  in  the  low-group  density  also. 

The  total  positive  ion  density  is  shown  in  Figure  3.9(d).  Since 
the  immobile  ions  show  the  same  gross  features  as  the  low-group  electron 
density,  it  is  clear  that  the  large  off-axis  peak  must  be  interpreted  in  terms 
of  the  history  of  the  ionization  rates.  However,  the  weak  ridge  in  the  low- 
group  electron  density  does  not  appear  in  the  ion  density,  and  thus  may  be  a 
flow  feature.  The  number  of  particles  in  the  high  group  is  not  large  enough 

for  the  ridges  of  Figure  3.9(c)  to  show  up  in  the  ion  density  (Figure  3.9(d)). 

It  seems  likely  that  they  are  similar  in  nature  to  the  ridge  in  the  low-group 
electron  density.  Further  support  for  the  flow  explanation  is  provided  by 
Figures  3.10(a)  and  (b).  These  show  the  low-group  and  total  ion-densities 
for  the  1-group  calculation  at  p/pQ  =  .005,  in  which  there  is  a  large  current 
enhancement  compared  to  p/pQ  =  .008  (see  Figure  3.7).  Here  the  ridge 
structures  in  the  electron  density  are  more  prominent,  but  still  have  no  strong 
counterpart  in  the  ion  density. 
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The  difference  in  positive  ion  density  radial  profiles  at  1  ns  is 
shown  clearly  in  Figure  3.11(a).  At  this  time  the  ion  density  is  much  greater 
than  the  beam  particle  density,  so  these  profiles  are  very  close  to  the  electron 
density  in  the  charge-neutralized  state.  The  effect  on  the  z-component  of 
plasma  current  was  shown  earlier  (Figure  3.8).  Figure  3.11(b)  shows  how  the 
relatively  small  difference  in  ion  densities  on  axis  due  to  time  delay  before 
01  ns,  becomes  very  large  for  a  few  tenths  of  ns  before  coming  together  again. 

For  ease  of  detailed  comparison,  electron  density  radial  profiles 
for  the  1-group  model  and  for  the  low  and  high-energy  components  of  the  3- 
group  model  are  shown  in  Figures  3.12(a)  -  (c).  The  large  temporary  off-axis 
hump  in  the  low-group  electrons  is  very  clear  at  0.6  ns  in  Figure  3.12(b). 

The  very  much  broader  density  profiles  produced  by  the  3-group  model  are  also 
very  obvious,  and  persist  out  to  10  ns  from  the  pulse  head. 

A  summary  of  electron  and  ion  number  density  comparisons  on  axis  is 
given  in  Table  3-1  below.  The  first  two  entries  of  the  third  column  show  the 
magnitude  of  the  time-delay  effect,  while  subsequent  entries  show  the  large 
differences  seen  in  Figure  3.11(b).  The  first  two  columns  give  an  indication 
of  the  importance  of  electron  transport.  At  0.3  ns  for  the  3-group  model,  only 
47%  of  the  plasma  electrons  ever  produced  on  axis  remain  there. 

TABLE  3-1 

ELECTRON  TRANSPORT  EFFECTS 

Time  (Ne/N.)  (Ne/N.) 

(ns)  1-Group  3-Group 


Figure  3.11(a).  Positive  Ion  Radial  Profile 
at  1  ns  (p/p  =  .008). 
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Figure  3.11(b).  Positive  Ion  Density  on  Axis 


3.12(a).  Low  Energy  Group  Electron  Fig.  3.12(b).  Low  Energy  Group  Electron 

Density  Profiles,  1-Group  Density  Profiles,  3-Group 

Model  (p/pQ  =  .008).  Model  (p/p  =  .008). 
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Fig.  3.12(c).  High  Energy  Electron  Density 
Profiles,  3-Group  Model 


Electric  Field  Effects 


Both  the  Ez  spike  and  the  final  decay  of  Er  occur  significantly  earlier 
in  time  for  the  1-group  model  (Figures  3.13(a)  -  (d)),  and  the  peak  field  strengths 
are  lower.  The  large  negative  oscillation  in  Er  at  1  Bennett  radius  is  especially 
prominent  in  the  3-group  calculation,  but  appears  to  damp  out  satisfactorily. 

The  negative  spike  in  Er  in  caused  by  an  overshoot  of  the  outward-moving  electrons 
as  charge  neutralization  is  finally  achieved  near  the  axis.  The  significant 
overshoot  is  consistent  with  the  fact  that  the  electron  plasma  frequency  is  much 
greater  than  the  collision  frequency  at  ^  0.3  ns.  The  large  off-axis  hump  in 
ion  density  at  0.6  ns  (Figure  3.12(b))  seems  to  be  associated  with  the  decelera¬ 
tion  of  the  low-group  electrons  by  the  reversed  E  ,  and  by  their  final  cooling 
through  the  peak  of  the  ionization  cross  section.  The  final  cooling  occurs 
considerably  earlier  off-axis,  as  shown  in  Figure  3.14(a).  The  temperature 
behavior  of  the  1-group  model  is  also  shown  for  completeness  in  Figure  3.14(b). 

3. 2. 5.4  Hall  Current  Effects 

Introduction 


Hall  currents  may  be  comparable  to  ordinary  currents  when  the  Larmor 
frequency  eB/mc  is  comparable  to  or  greater  than  the  momentum  transfer  frequency. 
Tnis  criterion  is  a  strong  function  of  both  ambient  density  and  electron  energy. 
The  3-group  model  allows  the  possibility  of  accounting  for  the  energy  dependence 
in  a  strongly  non-Maxwell ian  plasma.  The  momentum  transfer  frequency  at  1  eV 
in  is  about  the  same  as  that  at  200  eV,  and  decreases  with  energy  beyond  50 
or  60  eV.  Thus,  if  a  substantial  fraction  of  the  electrons  have  energy  in  excess 
of  200  eV,  their  low  momentum  transfer  frequency  may  result  in  a  significant  Hall 
current,  out  of  proportion  to  their  fractional  number  density. 


*o 

(a )  3-Group  Model . 
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1-Group  Model 


Figure  3.14(a)  and  (b).  Electron  Temperature  (0-1  ns)  ./; 


=  .008. 
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Detailed  Discussion 


It  is  clear  from  Figure  3.6  discussed  in  Section  3.2.3  above  that  Hall 
currents  are  important  to  the  pinch  force.  In  the  calculation  at  p/p  =  .008, 
even  the  low-energy  electron  group  achieves  high  velocity  for  a  short  time, 
and  contributes  substantial  Hall  current.  However,  this  is  generally  masked 
by  the  E^-driven  return  current.  In  the  far  wings,  however,  and  at  early 
times  near  the  axis,  the  radial  outflow  driven  by  Er  is  turned  forward  by  the 
magnetic  field  and  the  net  plasma  current  density  is  forward.  The  amount  of 
current  involved,  however,  is  insignificant  in  magnitude  compared  to  the  beam 
current. 

When  the  high-energy  electron  group  from  the  3-group  model  is  examined 
separately,  it  shows  quite  large  effects.  The  z-component  of  plasma  current 
for  the  high  group  is  shown  at  various  times  in  Figure  3.15(a).  Beyond  about 
1.5  Bennett  radii  it  is  always  directed  forward  with  the  beam.  Nearer  the  axis, 
the  strong  Ez  field  keeps  the  current  going  backwards,  although  by  10  ns  there 
is  relatively  little  backward-moving  current  in  the  high  group.  A  small  part 
of  this  forward-going  current  is  due  to  the  fact  that  high-group  electrons 
produced  directly  by  the  beam  are  injected  with  forward  velocity.  However, 
most  of  it  is  due  to  the  Hall  force  acting  on  the  radial  out-flow. 

The  effect  of  Hall  current  in  the  high  energy  electron  group  is  shown 

in  Figure  3.15(b).  At  1  ns,  the  forward-goi ng  current  near  1  Bennett  radius 

2  2 

amount  to  almost  2  kA/cm  ,  compared  to  the  5  kA/cm  of  low  group  plasma 
current  going  backwards  at  the  same  radial  position.  The  result  is  a  very 
significant  reduction  of  the  total  return  current. 

It  is  important  to  note  that  at  late  times  :  >  0.5  ns,  the  radial  flow 
which  is  turned  forward  by  the  magnetic  field  is  driven  not  primarily  by  Er, 
but  by  the  radial  pressure  gradient  of  the  high-group  particles.  Yu  (Ref.  10) 
first  recognized  this  as  a  possible  irportant  component  of  the  total  current. 


Figure  3.15(a).  Current  Density  J  of  High  Energy 
Electrons  (p/p  =  .008). 


Figure  3.15(b).  Effect  of  High  Energy  Electrons 

on  Total  Plasma  Current  Density  J 
at  1  ns  (.  /,  -  .008) . 


From  the  equation  of  motion  it  is  clear  that  a  density  gradient  is  equivalent 
to  a  radial  electric  field  of  magnitude 


\ 
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E  -  -  T  ■cv1n  N  volts/cm. 
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For  T  -o  10  KeV  as  suggested  by  Yu  (and  used  for  the  high  group  in  these  cal¬ 
culations),  and  a  density  gradient  with  scale  ^  1  cm,  this  corresponds  to 
Er  %  10  KeV/cm.  Assuming  that  the  system  as  a  whole  approaches  a  steady  state 
solution  of  Maxwell's  equations,  Jr  must  -*■  0.  In  the  presence  of  the  Hall- 
force  terms,  the  radial  Er  adjusts  itself  to  shut  off  the  net  current  driven 
by  v  B  ,  E  ,  and  the  pressure  gradient.  In  the  present  context  it  seems  that 
zero  net  radial  current  could  perhaps  be  achieved  by  a  balance  between  an  in¬ 
ward  flow  of  very  low  energy  electrons  and  an  outward  flow  of  a  much  smaller 
number  of  fast  electrons.  Because  of  the  differences  in  vm  for  the  two  streams, 
the  corresponding  Hall  currents  would  not  cancel  exactly  in  this  case.  Under 
some  conditions  such  counterflows  may  be  limited  by  instabilities. 

3. 2. 5. 5  Delta  Ray  Effects  at  Low  Density 

Although  the  steady  state  5 -ray  current  is  independent  of  air  density, 
the  time  required  for  the  current  to  reach  its  maximum  value  depends  strongly 
on  the  density.  There  is  also  a  significant  dependence  on  the  beam  energy  if 
only  relativistic  particles  (-  1  MeV)  are  included,  as  shown  in  Figure  3.16. 

The  pulse  length  used  in  this  calculation  was  100  ns.  (A  very  steep  drop  in 
5 -current  at  1  atmosphere  on  a  timescale  of  10  ns  beginning  just  before  the 
pulse  ends  is  not  shown.)  It  is  clear  that  at  densities  as  low  as  .01  atmosphere 
the  5 -current  represents  a  very  small  increment  to  the  beam  current  over  the 
first  10  ns  considered  in  the  calculations  presented  above.  Before  0.1  ns, 
the  5 -current  is  comparable  to  the  plasma  current  in  the  z-direction,  but 
both  are  very  much  smaller  than  the  beam  current.  It  seems  unlikely  that  the 
relativistic  part  of  the  beam  secondary  cascade  can  be  important  at  early 
times  for  air  densities  as  low  as  .01  atomsphere  However,  for  propagation 
in  a  density  channel,  the  interaction  with  the  high-density  walls  may  be  very 
important. 
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3. 2. 5. 6  Comparison  with  Ohm's  Law 


Both  the  1-  and  3-group  models  use  equations  of  motion  to  describe 
the  electron  flows.  In  general,  the  difference  between  the  currents  calculated 
this  way,  and  by  using  an  Ohm's  law  relationship  involving  only  electron  density 
and  momentum  transfer  frequency,  is  large  at  early  times.  The  duration  of  the 
deviation  from  Ohm's  law  increases  as  the  gas  density  decreases.  Figures  3.17(a) 
and  (b)  show  the  actually-computed  plasma  current  compared  with  what  would  have 
been  obtained  from  an  ohmic  calculation  using  the  instantaneous  values  of  Ng 
and  vm.  The  results  are  shown  for  the  1-group  calculation  at  p/pQ  =  .008,  but 
comparable  effects  are  present  for  the  3-group  case.  At  0.1  ns  the  ohmic  current 
shows  an  extremely  strong  Hall  effect  due  to  the  large  radial  field;  the  current 
given  by  the  equation  of  motion  is  very  much  smaller,  and  the  Hall  effect  does 
not  show  on  the  linear  scale.  At  1  ns,  the  two  methods  of  calculating  the 
current  give  virtually  identical  profiles,  as  shown  in  Figure  3.17(b). 

3. 2. 5. 7  Concluding  Remarks 

The  calculations  discussed  above  show  the  magnitude  of  the  effects  to 
be  expected  when  low-density  phenomena  are  taken  into  account.  All  the  commonly 
used  assumptions  listed  in  Section  3.2.2  above  have  significant  impact  on  the 
results  at  p/p  =  .01.  (The  current  carried  by  relativistic  beam  secondaries 

is  important  at  low  densities  only  if  the  pulse  is  long  enough  for  the  full 
build-up  to  occur. ) 
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Figure  3.17(a)  and  (b).  Comparison  of  Ohmic  Calculation  of 

Plasma  Current  Density  with  Equation 
of  Motion  Calculation. 
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4.0 


NONLINEAR  FIELD  ALGORITHM  DEVELOPMENT 


4.1  Outline  of  the  MAHD  Approach 

The  MAHD  (for  "Moderate  Amplitude  Hose  Displacement")  e.m.  field 
algorithms  were  developed  for  application  to  a  phenomenological  model  of  elec¬ 
tron  beam  hose  dynamics  In  the  nonlinear  regime.  The  value  of  a  phenomeno¬ 
logical  hose  dynamics  model  In  the  linear  (small  hose  displacement)  regime  has 
been  well  demonstrated  by  codes  such  as  PHLAP  (Ref.  1);  the  extension  to  the 
nonlinear  regime  Is  critical  for  studying  non-linear  saturation  and  propaga¬ 
tion  In  non-uniform  ambient  density  and  conductivity  channels.  PHLAP  and 
other  similar  models  approximate  some  of  the  physical  details  Incorporated  In 
elaborate  particle  simulations,  but  are  fast-running  and  economical.  With 
careful  calibration  against  comprehensive  simulation  codes  (Ref.  7),  phenome¬ 
nological  models  make  parametric  studies  over  extended  ranges  of  beam 
characteristics  more  feasible. 

The  MAHD  electromagnetic  (e.m.)  field  algorithms  are  required  to 
be  fast  In  execution  and  sufficiently  accurate  for  the  applications  of 
Interest.  All  field  solution  methods  presently  In  use  (DYNASTY,  Ref.  14, 
DYNADISC,  Ref.  15,  and  IPROP,  Ref.  7)  employ  a  modal  expansion  of  the  beam's 
fields  about  some  axis  of  symmetry;  MAHD  Is  not  an  exception.  Because  an 
accurate  representation  of  the  fields  for  a  beam  that  Is  strongly  displaced 
(say,  by  more  than  a  Bennett  radius)  from  the  nominal  propagation  axis 
requires  an  extremely  large  number  of  modes,  validity  for  moderate  displace¬ 
ments  only  Is  claimed  for  the  MAHD  algorithms.  This  limitation  Is  shared.  In 
practical  terms,  by  other  schemes  using  the  modal  expansion  approach  as  well. 

The  present  study  considered  two  distinct  physical  approximations 
In  the  formulation  of  the  field  equations  .  A  first  approximation  considered 
was  the  so-called  "frozen  approximation".  A  further  approximation,  developed 
by  Lee  (Ref.  16),  neglects  the  effects  of  the  Ez~  and  Hz~terms  that  appear  In 
the  field  equations.  One  has  the  further  option  of  considering  equations 
written  In  terms  of  potentials  or  of  dealing  directly  with  the  electron. agnetlc 
fields.  In  this  study,  we  compare  three  algorithms  as  follows: 


-82- 


o  Potential  equations  In  Lee's  approximation, 
o  Field  equations  In  the  frozen  approximation,  and 

o  Field  equations  In  Lee's  approximation. 

The  basic  solution  technique  employed  In  the  MAHD  algorithms  Is  as 

follows:  In  all  formulations  of  the  field  (and  potential)  equations  for  the 

modal  components,  coupling  between  modes  occurs  In  the  field  equations  because 

of  the  angular  dependence  of  ^ .  For  example,  the  m-th  mode  of  the  conduction 

current,  { oE }  ,  Is  just  the  coefficient  of  cos(mO)  In  the  expansion  of  the 
m 

product  of  the  two  series  for  o  and  for  E.  Besides  contributions  from  the 
mode  of  Interest,  o  E  ,  other  terms  a.  E,  appear,  where  |k  +  l|  =  m.  This 
coupling  requires.  In  principle,  that  the  equations  for  the  amplitudes  for  all 
modes  be  solved  simultaneously.  We  have  here  Instead  employed  a  scheme  close 
to  that  used  by  Godfrey  In  IPROP  (Ref.  7).  In  this  scheme  the  modally- 
resolved  set  of  field  equations  are  solved  In  mode  by  mode  subsets.  Contribu¬ 
tions  to  the  current  from  cross  terms  other  than  of  are  regarded  as  approxl- 
mately  known  and  fixed  for  the  m-th  mode  calculation.  Cross-term  contribu¬ 
tions  are  refined  In  an  Iterative  procedure  Involving  repeated  sweeps  through 
solutions  for  sequential  mode  values.  The  Iterative  scheme  avoids  the  labor 
In  the  alternate  approach  of  solving  the  very  large  set  of  equations  generated 
by  solving  simultaneously  for  the  field  amplitudes  In  all  modes.  We  have 
found  that  only  a  few  passes  are  required  (two  to  three,  at  most)  for  good 
accuracy;  fast  computer  execution,  and  a  nearly  linear  Increase  In  execution 
tine  with  the  number  of  modes  carried,  result. 

The  remainder  of  this  section  will  provide  details  of  the  three 
algorithms  studied  and  of  the  Iterative  solution  scheme  that  they  all  used.  A 
standard  test  problem,  formulated  by  F.  Chambers  of  LLNL  and  G.  Joyce  of  NRL, 
was  used  for  comparison  of  the  three  algorithms  among  themselves  and  with 
other  solution  codes;  results  of  those  comparisons  will  also  be  given. 
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4.2  Formulation  of  the  MAHO  Algorithms 

The  MAHD  algorithms  arc  formulated  with  the  frozen  approximation 
and  the  modal  expansion  as  a  common  basis.  Given  fields  In  the  laboratory 
cylindrical  coordinates  (r'»  6 • ,  z*,  t'),  we  make  the  usual  Galilean  trans¬ 
formation  to  beam  coordinates 


* 


(4.1) 


where  the  coordinate  €  measures  "beam  time"*  or  time  elapsed  since  the  head  of 
the  beam  (at  5=0)  passes  an  observer  at  some  value  of  z'.  The  frozen 
approximation  follows  Immediately  from  a  transformation  of  the  field  equations 
to  beam  coordinates  and  the  assumption  that*  In  beam  coordinates,  fields  vary 
slowly  with  z;  that  Is: 

J_  «  _9_  (4.2) 

3z  95 


Specific  sets  of  field  equations  will  be  shown  later  In  the  section. 

The  modal  expansion  Is  simply  a  Fourler-serles  representation  of 
the  dependence  of  the  various  field  quantities  on  the  transverse  coordinate  6. 
We  have  assumed,  to  simplify  Initial  development  work,  that  the  (x,z)-plane  Is 
a  symmetry  plane  for  the  fields,  corresponding  physically  to  hose  displace¬ 
ments  In  the  x-dlrectlon  only,  and  that  the -  and  r-components  of  the  beam 
current  can  be  Ignored.  With  these  assumptions,  the  primary  beam  current 
J  (r»fc,5)  Is  given,  for  example,  by: 

N 

Jz(r,  ,-)  =  y  *  Jzm(r,r)  cos  it:  (4*3) 

m  =  0 


k«  lit  "  «  i'J 
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The  beam-induced  conductivity  o»  the  e.m.  potentials  A  and  4>»  and  the  e.m. 
field  quantities  Er#  Ez#  and  HQ  are  also  even  (symmetric)  functions  of  6,  and 
are  given  by  analogous  cosine-series  expansions.  The  remaining  field  quan¬ 
tities  Hr»  H z*  and  E  rj  are  odd  (antisymmetric)  functions  of  6,  and  are  given  by 
a  sine-series  expansion;  for  example: 


Hz(r,e,0  = 


Hzm(r,c)  sin  me 


(4.4) 


m  =  1 


All  three  algorithms  use  essentially  Identical  methods  to  resolve 
the  (r»6)  sources  Into  modes  —  as  generally  required  by  the  algorithms  —  and 
to  form  the  modal  amplitudes  associated  with  the  products  of  two  modal  expan¬ 
sions  (as  required  by  the  Iterative  scheme).  We  briefly  summarize  those 
methods  next. 

The  present  algorithms  all  use  straightforward  Fourier  series 
expansion  formulae.  It  was  felt  that#  for  small  numbers  of  modes#  the  speed 
advantage  of  FFT  techniques  was  probably  not  significant.  The  modal  resolu¬ 
tion  procedure  does  In  fact  typically  use  only  about  5%  of  the  total  time 
spent  by  the  solution  algorithms  as  presently  programmed  (see  Sec.  4.4). 

By  our  convention,  an  even  function  n(fi)  Is  expanded  as 


+  V'  o  cos  m6 

L*t  m 


(4.5) 


where 


°o  =  T  f 

J  0 


■(fi)do 


(4.6) 


cm  it 


o(e)  cos  me  do 


m  >  1 


(4.7) 


(The  Integrals  In  the  above  equations  are  evaluated  using 
Simpson's  (3-polnt)  rule  In  a  direct  numerical  Integration.  If  li  Is  the 
highest  mode  number  to  be  used,  2N+1  points  in  the  Interval  CO# tt 3  are 
typically  used  In  the  Integration.) 

All  of  the  MAHO  algorithms  require  re-expanding  the  products  of 
two  modal  series  Into  modes.  Assume  that  a(e)  and  b ( G )  are  even  functions 
representable  by  a  finite  cosine  series  (maximum  mode  number  N),  and  that  c(0) 
and  d(6)  are  odd  functions  given  by  sine  series  of  the  same  length  N.  We 
Ignore  modes  of  higher  order  than  N  In  the  product. 

Even-Even  Products.  The  product  a(e)b(6)  Is  also  even,  and  the 
cosine  series  coefficients  are  as  follows: 


K 


‘obo  +  2  S  aibi 


(4.8a) 


m  >.  1: 


i  IWT1 

=  aobm  +  ambo  +  7  L  aTV(+  7  Y,  +  am+.b.)  “».8b) 


The  first  summation  on  the  right  Is  dropped  If  m  <.  1,  and  the  second  Is 
dropped  If  m  =  N. 


Odd-Odd  Products.  The  product  c(e)d(e)  Is  even,  and  the  cosine 
series  coefficients  of  are  as  follows: 


m  =  0: 


(4.9a) 


m  >  1: 


U  - 

(  )m 


m-1 


N-m 


2  2  C£dm-£  + 


£  =  1 


2  S(c£dm+£+  cm+£d£} 


i  =1 


(4.9b) 


As  before,  the  first  summation  on  the  right  Is  dropped  If  m  <  1,  and  the 
second  Is  dropped  If  m  =  N. 

Even-Odd  Products.  The  product  a(e)c(e)  Is  odd,  and  the  sine- 
series  coefficients  for  mil  are  as  follows: 


W 

i  L 


a  c  +  i 
o  m  2 


m-  1 
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£  =1 
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Summations  ere  dropped  for  m  =  1  and  m  =  N  as  In  the  previous  cases. 

As  Indicated  earlier,  we  employ  Iteration  schemes  In  which  (for 

example)  conduction  current  In  a  particular  mode,  { a E }  ,  Is  separated  Into  the 

m 

principal  contribution  from  the  mode  of  Interest,  o  E  and  other  contrlbu- 

o  m 

tlons.  We  thus  define  the  function  Pcc(a,b,n)  as  representing  the  summed 


terms  In  the  m-th  mode  of  the  product  of  the  two  cosine  series  for  the  even 

functions  a(0)  and  b(e),  except  for  the  term  a  b  : 

r  o  m 


Pcc(a,b,m)  =  {ab}m  - 

Inspection  of  Eq.  (4.8)  above  shows  that  Pcc(a#b,m)  as  defined  may  contain  the 

hlqh-order  term  o_  b  ,  but  contains  no  other  low-order  terms  Involving  b  .  We 

also  define  the  analogous  function  P  (a,c,m)  to  separate  a  c  and  other  con- 

3  cs  r  o  m 

trlbutlons  to  the  m-th  mode  of  the  product  of  the  cosine  series  a ( 6 )  and  the 
sine  series  c(0).  Note  that  In  this  case#  a(e)c(e)  Is  an  odd  function  of  0» 
and  Is  expanded  In  a  sine  series;  there  Is  no  m=0  component  In  the  product 
expansion: 


Pcs(a,c,m) 


{acl  -  a  c  (m  >  1) 
m  o  m  —  ' 


(4.12) 


Finally,  simply  define  P  (c,d,m)  as  the  m-th  mode  amplitude  In 

ss 

the  (cosine)  series  expansion  of  the  product  of  the  odd  functions  c(G)  and 
d  ( f. ) . 


Pss(c,d,m)  =  {cd}^ 


(4.13) 


Explicit  expressions  for  P  ,  P  ,  and  P  follow  Immediately  from 
r  r  cc  cs  ss 

Eqs.  (4.8)  -  (4.10). 

4.2.1  Reduced-Potential  Algorithm  MAHD1 

The  first  algorithm  we  will  consider  Is  based  on  Lee's  simplifi¬ 
cation  (Ref.  16)  of  the  frozen  approximation  for  the  e.m.  potential  equations. 
The  assumptions  leading  to  the  final  set  of  equations  are  equivalent  to  neg¬ 
lecting  beam-time  derivatives  of  Ez  and  Hz  In  the  e.m.  field  equations.  It  is 
unnecessary  to  find  the  transverse  component  of  the  vector  potential,  A^.  It 
is  convenient  to  define  the  remaining  two  potential  components  of  Interest  as 
A  and  :,  where  :  is  the  usual  scalar  potential,  and 


88- 


A  =  A  -  d> 
z  Y 

where  Is  the  axial  component  of  the  usual  vector  potential. 

In  RMKS  units#  Lee's  equations  for  A  and  41  are 


si(A  +  <  >  =  vl'  Vz 


V1  H  =  ^lZoo)  *  +  Z0oVi4> 


(4.14) 


(4.15) 


(4.16) 


(A  and  4>  In  the  above  equations  are  expressed  In  units  of  volts;  the  beam-time 
coordinate  £  Is  expressed  In  meters.  ZQ  Is  the  free-space  Impedance  of 
376.7  ohms#  and  the  units  of  the  conductivity#  o#  are  mho/m.)  The  e.m.  fields 
are  found  from  the  potentials  via  the  following  relations: 


?A 


V 


Hz  =  7  x 


H  =  -  u  x  V  a, 
1  z  i  2 


(4.17) 


(The  units  of  E  and  H  In  the  RMKS  form  above  and  to  be  used  hereafter  are 
volts/meter;  the  more  usual  RMKS  H-fleld  differs  from  our  usage  by  a  factor  of 
ZQ.)  As  noted  earlier#  A^  and  Hz  are  not  calculated  In  this  formulation.) 

Equations  (4.15)  and  (4.16)  are  recast  Into  finite-difference  form 

by  resolving  them  Into  m-th  modes  (0  <  m  <  N).  Potential  quantities  are 

defined  on  a  set  of  radial  grid  points  r.  (1  <  j  <  ND)  and  at  the  beam-time 

J  K 

points  5.|  (1  =  1#  2#  ....).  The  differencing  Is  centered  at  (1  +  1/2, j)»  l.e., 
on  the  radial  grldpolnts#  but  midway  between  ^  and  An  arbitrary  radial 

grid  structure  (In  terms  of  spacing)  Is  assumed#  and  3-polnt  Lagrange  dif¬ 
ferentiation  formulae  are  used  generally  to  evaluate  first-  and  second-radial 

derivatives  at  r.. 

J 


. 

'  N 


-F.9- 


Equation  (4.15)  Is  readily  differenced.  The  lefthand  term  Is 
represented  as  follows:  for  the  m-th  mode. 


\  2  ( 
Vj(A+*){ 


?  2 
a  1  J_  m_ 

.2  r  3r  "  2 

9r  r 


(A  +  6  ) 
v  rri  Ym' 


(4.18) 


Radial  derivatives  centered  at  Tj  In  the  above  equation  are  simply  expressed 
In  terms  of  potentials  at  the  points  rj_^»  rj>  and  rj  +  ^* 

The  first  term  on  the  r.h.s.  of  Eq.  (4.15)  Is  separated  (as 
suggested  at  the  beginning  of  4.2)  Into  two  portions: 


Z  a  4^  +  P  (Zc,~.ni) 
0  0  3t  CC  0 


(4.19) 


The  quantity  z0a0OA  /H)  Is  represented  as  a  simple  difference,  centered  at 

(i+l/2,j).  The  function  Pcc  (a,b,m)  represents  the  summed  terms  In  the  m-th 

mode  of  the  product  of  the  cosine  series  for  the  functions  a  and  Jz,  except  for 

the  term  a  b  .  P  contains  no  other  low-order  b  mode  terms,  and  Is  con- 
o  m  cc  m 

sldered  a  known  quantity  In  the  solution  of  the  m-th  mode  difference  equations 
derived  here.  Specific  expressions  for  Pcc  and  other  related  product  expan¬ 
sions  are  found  from  Eqs.  (4.8)  -  (4.13). 

Collecting  the  m-th  mode  fields  on  the  left,  Eq.  (4.15)  finally 

yields 


r>  pin 

d  ( A  +<{,)-  Z  C  ;  •— 
m  m  m  o  o  V" 


P(Znc,  4r  ,  m)  -  Z  J 
cc  o  d.r,  o  zr 


(4.20) 


where 
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(4.21) 
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1  J_  nr 
r  3r  " 


The  finite-difference  representation  of  Eq.  (4.16)  Is  obtained  In 
an  analogous  manner  to  Eq.  (4.15).  The  lefthand  term  of  Eq.  (4.16)  Is  written 
•  (see  Eq.  (4.18)) 


3A  ( 

)  '-1-  \  m 


(4.22) 


with  finite-difference  repesentatlons  of  the  derivatives  centered  as  explained 
above. 

The  first  term  on  the  r.h.s.  of  Eq.  (4.16)  Is  a  vector  product: 

(z  c)  •  v  ($)  =  z  i£  ii  +  I  .  I  (4.23) 

0  1  o  ar  3r  r  aS  r  36 

The  radial  derivative  product  Is  a  product  of  cosine  series#  and  Is  handled  as 
before: 
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3c  os  / 

3r  3r  k 
)m 


=  Z  o 
0  0 


+  p 

cc 


3i)> 

ar 


m, 


(4.24) 


The  e-derlvatlves  In  Eq.  (4.23)  generate  a  Fourier  sine  series  for  each  factor 
In  their  product;  for  example 


n.  =  1 


ma  sin  rr •• 
m 


(4.25) 


Because  the  sine  series  does  not  Include  a  monopole  term*  b  does  not  contrl- 

m 

bute  In  low  order  to  the  m-th  mode  of  the  sine-series  product. 

The  rightmost  term  In  Eq.  (4.16)  Is  handled  as  the  previous 
cosine-series  products: 


I 2  ,v2  I 

,w, 


p  p 

=  Z  cj  dif.  +  P  (Z  o  ,  V,d>.  m) 
o  o  nTm  cc  '  o  o’  j/ *  ' 


(4.26) 


Collecting  the  m-th  mode  fields  on  the  left,  Eq.  (4.16)  finally  yields 


.2  _ m 

m  35 


0  3o  36 

7  7  on 

Zo°o  ”  Zo  3r  3r 


P  (Z  o,  ?  ,  m)  +  P  (  —  ,  m 

cc'  o  3r  ss  \  r  36  r  30  ’ 


(4.27) 


+  Pcc  (Zo0’  Y’  m) 


The  finite-difference  forms  of  Eqs.  (4,20)  and  (4.27)  are  applied 
to  each  of  the  Interior  grldpolnts  Tj  (2  <  j  <  NR_^).  Two  Inner  and  two  outer 
boundary  conditions  must  be  Imposed  to  complete  the  set  of  equations  for  the 

NR  unknown  values  cf  A  and  We  adopt  the  standard  boundary  conditions. 


at  r  =  0: 


3A  (o  96  (o) 
_ o =  n  — -2 

3r  ’  3r 


Am(o)  =  0,  6n(o)  =  0,  for  m>l; 


(4.28) 


atr=R  :  A  (R  )  =  0;  t_(R  )  =  0,  for  all  m 
max  mv  max'  ’  mv  max'  ’ 
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The  band  matrix  representing  the  set  of  2  NR  equations  for  the  A's 
and  <Ts  has  a  bl-trldlagonal  structure;  the  equations  can  be  solved  readily 
via  any  of  several  schemes,  Including  straightforward  Gaussian  elimination. 

The  computational  cycle  for  a  particular  field  tlmestep  Is 
Illustrated  In  Fig.  4.1,  which  outlines  the  major  aspects  of  the  Iterative 
approach.  Equations  (4.20)  and  (4.27)  (with  the  boundary  conditions  of 
Eq.  (4.28))  are  solved  for  potential  amplitudes  at  successive  m-values.  At  a 
given  point  In  a  pass  through  the  m-values,  the  arrays  from  which  the  P  and 

vv, 

Pss  terms  are  calculated  Involve  recently-updated  mode  amplitudes  for  m-values 
less  than  the  mode  of  current  Interest,  and  "last-pass"  amplitudes  for  modes 
equal  to  and  higher  than  the  current  mode.  Updated  amplitudes  are  entered  for 
use  in  the  Pcc  and  Pgs  computations  as  they  are  generated.  Only  a  few  Itera¬ 
tions  (two  to  three  passes  at  most,  for  the  test  case  considered)  appear  to  be 
required  for  accurate  results.  (A  discussion  of  comparisons  will  be  given  In 


Sec .  4.4). 


4.2.2 


Frozen-field  Algorithm  MAHD2 


The  next  algorithm  to  be  considered  Is  a  calculation  that  treats 
e.m.  fields  directly  (rather  than  through  potential  functions)  and  with  no 
simplifications  beyond  the  frozen  approximation.  Inserting  the  transforma¬ 
tions  (4.1)  and  the  frozen  approximation  (4.2)  Into  Maxwell's  equations 
Immediately  yields  the  basic  field  equations  treated  In  MAHD2  and  (with 
further  approximations)  MAHD3: 


:-E  ,  >H 

uV(ViEz  ♦!?*>*'  — 


vi-i-",)  =  Vz4  — 
-Hi 

H;  '  —  >  '  Vi4 


(4.29) 


v**A 


Figure  4.1.  Iterative  computation  cycl 


The  RMKS  Maxwell  equations  are  readily  decomposed  Into  modal 
equations;  maintaining  the  symmetry  assumptions  outlined  at  the  beginning  of 
the  Section  (even  angular  symmetry  for  o»  J  #  E  »  E  .  and  H  ;  qi id.  angular 

Z  Z  r  u 


symmetry  for  H^,  Hr#  and  E0)#  we  find: 


-  rr  (  rE  )  +  -  E 
r  3r  r  rm 


(4.30a) 


(4.30b) 


-  —  E 


r  zm 


(4.30c) 


I  3  tru  \  m  u  _  °  zw  7  (  ) 

r  sr  ^vm  "  r  Hrm  ~  3£  Zo  |Jz  +  °^z  j 


(4.30d) 


V +  zo  K 


(4.30e) 


+  —  H_ 


—  +  Z  <  n  F  ( 


(4.30f ) 


Note  that  the  mode  amplitude  { aE0 } m  In  Eq.  (4.30e)  Is  for  an  odd-symmetry 
mode;  the  adjacent  equations  are  for  even-symmetry  modes. 


C 


Because  we  have  not  Included  any  components  of  the  beam  current 

but  J  ,  the  monopole-case  (m=0)  fields  collapse  to  the  TM  set  E  ,  E  .  and  H„. 
z  z  r  o 

Solutions  for  this  axlsymmetrlc  case  are  quite  standard,  and  will  not  be 
considered  In  detail  here.  We  will  concentrate  on  the  (m>l)-modes  in  the 
following  discussion. 

Equations  (4.30a-e)  are  differenced  In  a  straightforward  fashion. 
First,  note  that  Eqs.  (4.30c)  and  (4.30f)  are  local  equations  In  the  sense 
that  they  do  not  Involve  spatial  derivatives.  Finite-difference  representa¬ 
tions  of  the  "c"  and  wf"  equations  are  thus  centered  at  ( 1+1/2, J)  and  written 
for  each  of  the  NR  radial  grldpolnts.  For  one  boundary  condition  treatment, 
special  forms  of  "c"  and  "f"  are  required  at  r^=0,  and  will  be  discussed  In 
connection  with  boundary  conditions  on  the  fields,  below. 

Equation  (4.30c),  differenced  with  the  Indicated  centering,  leads 

to  a  relation  for  E  (1+1, J)  In  terms  of  Eft  (1+1, j)  and  H  (1+1, J): 

zm  om  rm 


Ez,(i 


+  l,j)  =  c 
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E,  (i  +  l,j)  + 


■m 


H  (i 

rm 


+  l.j) 


'2j 


(4.31) 


where 


■2i 


*  c 


Ij 


(4.32) 


Equation  ( 4 . 30 f )  Is  formally  Integrated  In  time  to  obtain  a 

relation  for  H  (1+1, J)  In  terms  of  E  (1+1, j)  and  H  (1+1,1): 
zm  rm  m 
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(4.33) 


W1  +  2,j) 


aErm{i*j) 
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+  n 

3t  r  z 


ZoPcc(o>Er>  m) 


yielding 


where 


HZn, ( 1  +  l.j)  =  F,  E„(1  ♦  l.j)  -  F2Hen(i  +  l.j)  +  F 


1  rm 


2  en' 


(4.34) 


a.  =  exp(-Z0o0A£) 


z0oA£ 


Flj 


F2j 


2r. 
_ J_ 

n^j 
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(4.35) 


F3j  =  F2j(H  (i>j)  +  Pcc(Zo0^,Er»m))  "  3jFljEr(iJ)  "  Hz(iJ) 


The  remaining  equations*  which  Involve  radial  derivatives*  are 
cast  into  finite-difference  forms  centered  at  ( 1+l/2»J+l/2) .  (Again*  special 
forms  of  Eqs.  (4.30a*  b»  d,  and  e)  are  written  at  (l+l/2*3/2)»  the  first 
radial  differencing  cell.  In  connection  with  field  boundary  conditions;  they 
will  be  discussed  later.)  Equations  (4.30a)  and  (4.30b)  are  differenced  In 
analogy  with  Eq.  (4.30c)  (see  Eq.  (4.31)  above),  while  Eqs.  (4.30d)  and 
(4.30e)  are  formally  Integrated  In  analogy  to  Eq.  (4.30f).  The  differenced 
form  of  Eq.  (4.30a)  (dropping  the  subscript  "m"  and  the  Index  n1+lw)  Is  thus: 

AnE  (j  +  1)  -  A2E  (j)  +  A3(Er(j  +  1)  +  Er(j))  +  H,(j  -  1}  +  H^j)  (4.36) 

-  r  h  s , 


where 


A  =^iJ± 
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fl  r  _ - 

h2  r.  r 
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A  = 
m3  r . 

J+’2 


(4.37) 


rhsA  =  Hz^»j  +  1)  +  "  AiEr.  (i  »J+1)  + 


-  A3(Er(i,j)  +Er(i,j+1)) 


Equation  (4.30b)  Is  differenced  analogously  to  (4.30a) ,  resulting  In  an 
equation  of  the  form: 


Er(j  +  1)  +  Er^j)  "  hcAJ 


(j  +  1)  -  Mj)  +  B,(Ez(j  +  1)  -  Ez(j)) 


(4.38) 


=  rhsr 


The  differenced  form  of  Eq.  (4.30d)  Is  obtained  by  a  formal 


Integration 


E  (i  +  1 ,  j  +4)  =  aE  (i ,3  +  j) 


(4.39) 


+  (1^1  I  -4~  (rH  )  -  —  H  -  Z  J  -P  (Z  a,  E  ,  m) 
Z  o  r  3r  ••  r  r  o  zn;  cs  0  z 
o  o 


•here  the  quantities  enclosed  In  square  brackets  are  evaluated  In  ordinary 
difference  form#  centered  by  ( 1+1/2# J+l/2) .  Expansion  of  this  equation  leads 
to 


Ez(J  +  D  +  E2(j)  '  B( j+-,i)A1H  (j+1 )  +  f.(j+la)A2H  (j) 
+  f  (j+:2)A3(Hr(j+l)  +  Hr  ( j ) )  =  rhSp 


(4.40) 


where 


rhsD  =  a(j+Js)(Ez(i,j+l)  +  Ez(i,j))  +  BU+^AjH  Ji  ,j  +  l) 

-  e(j+,i)A2H  (i,j)  -  B(j+*s)A3(Hr(i,j+l)  +  Hr(i,j))  (4.41) 

-  2i5(j+Js)(Z0ACJ2  +  PccCZoO.Ez#"1)) 


A  precisely  analogous  procedure  Is  applied  to  Eq. 
formal  Integration: 


(4.30c), 


beginning  with  a 


a-*) 

Z  c 


P 

cs 


(4.42) 


m) 


We  finally  obtain  from  (4.42)  and  equation  of  a  form  analogous  to  (4.38)# 

Involving  H  ,  E.  #  and  H  . 

3  r  r  z 


Equations  (4.31)  and  (4.34)  are  used  to  eliminate  Ez  and  Hz  from 

(4.36)  -  (4.39)#  generating  In  principle  an  Incomplete  set  of  4(Np-l) 

equations  In  the  4N^  unknown  fields  Erm(i  +  l#j)»  E^H  +  l#,)).  1  + 1  #  J )  #  and 

H  (1+1,1).  Two  outer  boundary  conditions  and  two  Inner  boundary  conditions 
•m 

are  required  to  complete  the  set  of  equations  and  make  It  soluble. 


i 


Simple  conducting-boundary  conditions  are  always  Imposed  at  the 

outer  wall  at  r..c  =  R  .  namely 
nk  max 


E  (R  )  =  0 

max7 


H  (R  )  =  0 

r  max7 


(4.43) 


At  the  Inner  boundary,  several  sets  of  boundary  conditions  have 
been  tried.  A  satisfactory  set  appears  to  be 


E  -  E  = 


0 

H  +  H_  =  0 
0 
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E.. 


H_  = 


(4.44) 


with  the  latter  two  equations  replacing  nc"  and  wf".  As  an  alternative  pro¬ 
cedure,  a  more  complicated  procedure  can  be  followed  at  the  Inner  boundary, 
r^  =  0.  Because  of  the  1/r-factors  In  Eqs.  (4.30c)  and  (4.30f),  It  Is 
necessary  to  look  more  closely  at  the  behavior  of  the  fields  as  r  approaches 

zero.  We  find  that,  near  r  =  0,  E  ^  rm,  H  rm;  and  the  transverse  fields 

zm  .  zm 

E  »  E_  ,  H  ,  and  H  all  vary  as  r™-  .  Within  the  first  radial  grid  cell, 
rm  rm  rm  rm  3 

we  thus  assume 


Ern(r-«> 


(4.45) 


and  so  on,  for  the  remaining  field  quantities.  We  find  that,  at  r  =  0,  the 
field  equations  (4  30c)  and  (4.30f)  are  replaced  by  equivalent  relations 
between  the  e's  and  h’s: 
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(4.46) 


3e 

~ 


35 


3e. 


+  Z  c  e 
o  o  r 


3h 

3 


—  +  mh_ 
r  z 


Vr 
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and  that  the  remaining  field  equations  yield  the  two  Independent  boundary 
conditions 


e  +  er  =  0 

h..  +  hr  =  0 


(4.48) 


The  remaining  equations  (4.30a,  b,  d,  and  e)  are  written  for  the  Interval 

r^  1  r  <  r^  In  terms  of  the  e’s  and  h*s>  but  are  otherwise  expressed  In 

finite-difference  form  with  the  handling  and  centering  described  above  for  the 

"usual"  situation  of  the  differencing  cells  between  r_  and  R  .  We  thus 

z  max 

obtain  a  set  of  4NR  equations  for  the  4(NR-1)  fields  at  rr  =  r^,  r^,  ...  r^R 

and  e  ,  e.,  h  ,  and  defined  at  r. .  The  axial  fields  E  ,  H  ,  e  ,  and  h  are 
r  6  r  e  1  z  z  z  z 

und  by  back-substitution  In  Eqs.  (4.31),  (4.34),  and  their  analogues, 
obtained  from  Eq.  (4.46)  and  (4.47),  once  the  solution  for  the  transverse 
fields  Is  complete. 

The  MAHD2  calcul atlonal  cycle  Is  carried  out  exactly  as  the  MAHD1 
scheme  discussed  In  connection  with  Fig,  4,1.  Arrays  used  to  determine  the 
conduction  current  appearing  In  the  various  modes  are  assumed  known  (except, 
of  course,  for  the  ooEm  contribution  separated  from  the  remainder)  and  updated 
as  the  sequential  solutions  for  mode  amplitudes  proceed.  The  solution  to  the 
monopole-mode  (m=0)  equations  simplifies  to  the  usual  axlsymmetrlc  field  equa¬ 
tion  solution,  and  Is  carried  out  In  a  separate  algorithm.  The  monopole  solu¬ 
tion  Is  Integrated  Into  the  iteration  procedure,  however;  the  contributions  to 

the  monopole  conduction  currents  {aE  )  and  {cE  }  are  evaluated  from  the  yet- 

r  c  z  o 

to-te  corrected  higher-order  mode  flelcs,  and  solutions  for  m  >  1  follow  the 
monopole  solution  in  the  iterative  cycle,  as  they  co  in  the  MAHD1  calculation. 


4.2.3  Simplified  Field  Algorithm  MAHD3  i 

( 

i 

In  order  to  extend  the  range  of  comparisons  possible,  a  variation  | 

of  MAHD2  was  also  constructed.  The  physical  assumptions  that  led  to  the 

further  simplifications  of  the  frozen  approximation  used  In  the  Lee  potential 

formulation  algorithm,  MAHD1,  amount.  In  field  language,  to  neglect  of  the 

terms  5H  /3£  and  3E  /3£that  appear  In  Eqs.  (4.30a)  and  (4.30d)  above.  ! 

z  z 

The  finite  difference  equations  solved  In  the  MAHD3  algorithm  were 
derived  by  modifying  the  Eqs.  (4.30a)  and  ( 4 . 3 Od )  as  noted  above.  In  pro¬ 
gramming  terms,  the  modifications  are  quite  small.  The  differenced  form  of  | 

Eq.  (4.30a)  shown  In  Eq.  (4.36)  Is  replaced  by: 

A-j  E  ( j  +  1)  +  A2E  (j)  +  A3(Er(j  +  1)  +  Er(j) )  =  rhs'A  (4’49) 


Since  no  longer  occurs  In  the  above  equation.  It  Is  used  directly  In  the 
equation  set  to  be  solved. 

Equation  (4.30d)  Is  no  longer  formally  Integrated,  and  Eq.  (4.39) 
Is  replaced  by: 

1  (4.50) 

J  ■  - 

i  ,  ,  * 

1  ^  O  *  v.  'o 


Boundary  conditions  and  the  handling  of  the  ( r=0)-vers1ons  of  Eq.  (4.30c)  and 
(4.30f)  (shown  In  Eqs.  (4.46)  and  (4.47),  respectively)  are  Identical  to 
MAHD2. 


4.3 


Auxiliary  Models  and  Methods 


The  algorithms  discussed  In  Sec.  4.2  above  are  Intended  to  Inter¬ 
face  In  a  reasonably  tidy  way  with  other  portions  of  a  complete  hose  dynamics 
model.  The  principal  Inputs  to  the  field  algorithm  are  thus  quite  simply  the 
two-dimensional  field  source  arrays  of  primary  beam  current  density  and  con¬ 
ductivity  that  reflect  the  spatial  <r,e)  variation  of  those  quantities  at  the 
currently  Interesting  5-step.  Such  field  source  arrays  may  be  defined  in  any 
way  that  suits  the  physics  of  the  desired  model  (e.g.»  analytic  formulae# 
particle  simulation  results*  or  phenomenological  models).  The  field  sources 
are  then  passed  on  to  the  field-advancement  scheme*  which  operates  Indepen¬ 
dently  of  the  details  of  their  origin.  At  present  only  simple  analytic  pre¬ 
scriptions  for  Jz  and  o  have  been  used  to  exercise  the  algorithms.  Details  of 
those  prescriptions  for  a  standard  test  problem  will  be  considered  In  the 
following  section. 

4.4  Test  Calculations  and  Comparisons 

The  algorithms  described  above  were  programmed  In  detail  for 
actual  machine  execution.  Since  Intercomparison  of  the  algorithms  (yia  a  yis 
their  accuracy,  speed,  the  effects  of  further  simplification  of  the  frozen 
approximation,  and  differences  between  field-  and  potential-  formulations  of 
the  e.m.  equations)  was  of  most  Interest,  a  simple  prescribed  model  of  beam 
current  and  conductivity  was  used  In  exercising  the  MAHD  algorithms. 

4.4.1  NRL  Standard  Nonlinear  Test  Problem 

The  beam  model  used  In  the  test  calculations  was  a  prescription 
developed  by  F.  Chambers  of  LLNL  and  G.  Joyce  of  NRL  as  a  standard  test 
problem  for  nonlinear  field  algorithms.  The  beam  current  density  distribution 
Is  assumed  to  have  a  Bennett  profile  characterized  by  af;  -Independent  Bennett 
radius  parameter,  r^.  The  total  primary  beam  current  1^(5)  Is  given  by  the 
analytic  form 
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Specific  parameters  for  the  Standard  Test  Case  (STC)  are  summarized  In 

Table  4-1. 

m 

4.4.2  Results  of  Calculations 

The  Standard  Test  Case  (STC)  problem  was  run  for  each  of  the  three 
9  algorithms  described  above.  Identical  radial  grlddlng  schemes  (60  grldpolnts 

on  an  expanding  mesh  to  Rwa-|-|  =  10  cm)  and  tlmestep  Intervals  (A£  =2x10  **sec) 
were  used.  Five  radial  modes  (m  =  0#  1,  2#  3#  4)  were  used#  and  three  Itera¬ 
tions  per  tlmestep  (probably  one  more  than  necessary*  see  below)  carried  out. 
C  In  calculations  run  to  5  =  1  nsec#  typical  elapsed  CDC-7600  cpu  times  were 

between  about  4  and  6  seconds#  for  all  of  the  algorithms.  The  MAHD1  algorithm 
was  also  reprogrammed  to  take  some  advantage  of  the  vectorlzatlon  capabilities 
of  the  Cray-1.  Calculations  to  %  -  2  nsec  with  MAHD1  on  a  Cray-1  used  a  total 
v  of  3.2  sec,  with  about  1  sec  used  In  the  output  routine.  Of  the  time  not  con¬ 

sumed  In  output  operations#  most  of  the  time  was  used  In  solving  the  120-row 
band  matrix  (240-row  for  MAHD2  and  MAHD3)  representing  the  differenced  equa¬ 
tions  for  the  unknown  potentlals/f lelds.  The  algorithms  appear  quite  fast  In 


TABLE  4-1 


SUMMARY  OF  STANDARD  TEST  CASE  PARAMETERS 
Beam  Current 


*b-pk  =  10' 000  amperes 
=  1  nsec 


r 

Ur 


30  cm) 


Beam  Current  Density 


B 


1  cm 


Beam  Displacement 


X 

0 

= 

1 

cm 

To 

= 

1 

nsec 

(<o 

30 

cm) 

T 

w 

s 

.25 

nsec 

% 

7.5 

cm) 

Conductivity 


bkg 


=  1.113  x  10-2  mho/m  (108  sec'1) 
=  350  mho/amp 
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comparison  with  the  time  (over  two  minutes)  required  for  the  same  problem* 
using  the  full-mode-set-ln version  approach. 


Elapsed  time  for  the  MAHD  algorithms  should  scale  with  the  product 
of  separate  linear  dependence  with  time  and  space  grid  size,  number  of  modes* 
and  number  of  Iterations.  In  the  STC  problem,  differences  between  I  and  2 
Iterations  led  typically  to  differences  of  about  .5%  In  the  calculatlonal 
results  at  1  nsec;  differences  between  2  and  3  Iterations  were  less  than 
.05%. 


Detailed  comparison  of  the  results  of  the  three  calculations 
showed  nearly  Identical  field  profiles  for  MAHD1  and  MAHD3.  Individual  field 
values  compared  to  better  than  2%  everywhere,  with  comparison  to  better  than 
.5%  for  most  fields  and  locations.  This  agreement  Is  encouraging  from  a 
programming  standpoint,  since  the  MAHD1  and  MAHD3  algorithms,  although  con¬ 
taining  the  same  (or  equivalent)  physical  assumptions*  are  quite  different  In 
analytic  and  programming  details.  Comparisons  with  a  few  data  points  avail¬ 
able  from  DYNASTY  results  for  the  STC,  as  quoted  In  Ref.  15  »  are  also  en¬ 
couraging.  The  results  concern  f lei d/potential  profiles  along  the  beam- 
displacement  (X-)axls  at  K  -  1  nsec,  the  time  of  maximum  beam  displacement  In 
the  STC  (see  Table  4-1).  Values  at  r  =  0  and  maximum/minimum  values  along  the 
outgoing  radius  at  6  =  0  are  tabulated  for  comparison  In  Table  4-2.  Dif¬ 
ferences,  which  are  not  large,  may  probably  be  attributed  to  differences  In 
time-  and  space-grlddlng  and  details  of  the  source  (current  and  conductivity) 
algorithms  In  DYNASTY  and  the  MAHD  codes.  More  detailed  comparisons  with  the 
DYNASTY-1  Ike  field  solver  In  the  DYNADISC  code  are  anticipated. 

Detailed  MAHD2/MAHD3  comparisons  are  shown  In  Figs.  4.2  through 
4.6.  The  figures  show  STC  results  at  £  =  .5  nsec  and  5=1  nsec  for  E  -  and 
H0-f1elds.  Monopole  and  dipole  modes  are  plotted  separately,  together  with 
the  "cumulative”  values  representing  actual  radial  profiles  along  the  direc¬ 
tion  of  beam  displacement,  6=0. 

Radial  E-flelds  are  compared  In  Fig.  4.2  at  K  =  0.5  nsec.  The 
monopole  fields  are  essentially  Identical  from  MAHD2  and  MAHD3,  but  dif¬ 
ferences  In  the  dipole  modes  are  fairly  large  and  significant  Inside  the  beam. 
A  similar  comparison  at  i  =  1  nsec  Is  shown  In  Figs.  4,3  and  4.4.  Sizeable 
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FIELD- IV/M) 


R(M) 


Figure  4.2.  STC  Er  profiles  at  •;  =  .5  nsec.  The  m=0  fields  for 

both  frozen  and  Lee  approximation  are  almost  identical 
Differences  in  the  m=l  modes  do  affect  the  cumulative- 
field  comparison  inside  the  beam,  however. 
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10 

RIM) 


Figure  4.3. 


Comparison  of  monopole  and  dipole 
Er  amplitudes  vs„  radius  at  t,  =  1 
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ampl itudes. 
nsec. 
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Figure  4.4.  Comparison  of  cumulative  fields  (-:=0)  at  •'  =  1  nsec. 

(Monopole  and  dipole  components  were  compared  in 
Fig.  4.3.) 
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Rim 

Figure  4.5.  STC  H. -profile  at  ■;  =  .5  nsec.  Monopole  components 
dominate  the  fields,  and  differences  in  the  m=l 
fields  have  little  effect  on  the  cumulative  profiles. 


differences  appear  there,  both  In  the  monopole  and  dipole  modes,  and  therefore 
again  In  the  cumulative  fields  (plotted  separately  for  clarity  In  Fig.  4.4)  as 
well.  Outside  of  the  beam  (here,  beyond  about  2  cm),  the  fields  again  are 
similar. 

Comparisons  for  the  magnetic  fields  are  typically  closer. 

Figure  4.5  shows  (-^-profiles  for  the  STC  at  C  =  .5  nsec.  The  relatively  small 
differences  In  the  dipole-mode  fields  do  not  significantly  affect  the  cumula¬ 
tive  fields.  At  C  =  1  nsec,  there  Is  essentially  no  difference  between  the 
MAHD2  and  MAHD3  results;  only  the  results  for  MAHD2  are  shown  In  Fig.  4.6, 
which  shows  the  relative  amplitude  for  all  five  modes  carried,  as  well  as  the 
summed  profile. 

Since  the  magnetic  fields  dominate  the  particle  dynamics  of  the 
beam,  the  E-fleld  differences  shown  In  the  calculatlonal  results  may  not  be 
critical.  On  the  other  hand,  breakdown  effects  In  the  beam  (not  modelled  In 
the  STC)  may  differ.  Further,  the  higher-mode  differences  are  more  Important 
for  dynamics  effects  than  the  (m=0)-f lelds,  so  that  further  comparison  cal¬ 
culations  would  be  of  Interest. 

4.5  Summary 

Three  alternative  formulations  of  the  nonlinear  e.m.  field 
solution  were  developed  and  programmed.  Two  formulations,  one  based  on  Lee's 
approximate  potential  equations,  and  another  using  an  equivalent  treatment 
explicitly  considering  e.m.  fields,  compared  extremely  closely  with  one 
another,  and  quite  well  with  a  few  data  points  available  from  a  DYNASTY 
calculation.  The  third  calculation  contained  the  unsimplified  frozen-field 
assumption  as  a  physical  basis  for  the  solution  formulation.  Moderate 
differences  In  the  fields  found  with  this  model  and  with  the  previous  two 
appear  to  exist.  All  of  the  models  use  Iteration  to  find  a  satisfactory 
solution  to  the  coupled-mode  aspects  of  the  problem;  the  Iteration  appears 
accurate  and  extremely  fast. 
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